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ABSTRACT 

This report deals with the development of a general mathematical model and solution 
methodologies for analyzing structural response of thin, metallic shell-like structures under 
dynamic and/or static thermomechanical loads. In the mathematical model, geometric as well as 
material-type of nonlinearities are considered. Traditional as well as novel approaches are reported 
and detailed applications are presented in the appendices. The emphasis for the mathematical 
model, the related solution schemes, and the applications, is on thermal viscoelastic and 
viscoplastic phenomena, which can predict creep and ratchetting. 

1. INTRODUCTION 

The prediction of inelastic behavior of metallic materials at elevated temperatures has 
increased in importance in recent years. The operating conditions within the hot section of a rocket 
motor or a modern gas turbine engine present an extremely harsh thermo-mechanical environment. 
Large thermal transients are induced each time the engine is started or shut down. Additional 
thermal transients from an elevated ambient occur, whenever the engine power level is adjusted to 


meet flight requirements. The structural elements employed to construct such hot sections, as well 
as any engine component located therein, must be capable of withstanding such extreme 
conditions. Failure of a component would, due to the critical nature of the hot section, lead to an 
immediate and catastrophic loss in power and thus cannot be tolerated. Consequently, assuring 
satisfactory long term performance for such components is a major concern for the designer. 

Traditionally, this requirement for long term durability has been a more significant concern 
for gas turbine engines rather than rocket motors. However, with the advent of reusable space 
vehicles, such as the Space Shuttle, the requirement to accurately predict future performance 
following repeated elevated temperature operations must now be extended to include the more 
extreme rocket motor application. These blades operate in severe thermal transients that result in 
large inelastic strains, and several types of behavior must be considered.The elevated temperatures 
can lead to thermal buckling and, in addition, creep can be expected to occur. Thus, a combination 
of thermal-creep buckling behavior leading to large deflections can be anticipated. Because of the 
cyclic character of the mechanical and thermal loads, progressive growth or ratchetting effects must 
also be considered. Thus, geometric and material nonlinearities (of high orders) can be anticipated 
and must be considered in the mathematical model. 

Consequently, the industry is concerned with the behavior of thin shell-like structural 
elements subjected to severe time-dependent thermomechanical loading. Such thin elements, 
including beams, rings, arches, plates and shells, are presenting generic types of components, 
which might be located within or adjacent to the hot section of a rocket or a gas turbine engine. 

The experience in the gas turbine engine industry indicates, however, that existing analytic 
tools are not sufficiently reliable to accomplish this task. State of the art methods for predicting hot 
section component behavior are generally not sufficiently accurate to perform extended use 
evaluations. 

Under this kind of severe loading conditions, the structural behavior is highly 
nonlinear due to the combined action of geometrical and physical nonlinearities. On one side, finite 
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deformation in a stressed structure introduces nonlinear geometric effects. On the other side, 
physical nonlinearities arise even in small strain regimes, whereby inelastic phenomena play a 
particularly important role. From a theoretical standpoint, nonlinear constitutive equations should 
be applied only in connection with nonlinear transformation measures (implying both deformation 
and rotations). However, in almost all of the works in this area, the two identified sources of 
nonlinearities are always separated. This separation yields, at one end of the spectrum, problems 
of large kinematics, while at the other end, problems of viscous and/or non-isothermal behavior in 
the presence of small strain. 

Because of the nature of the causes, special care is needed in the selection or development 
of a constitutive law that includes time and temperature effects. Although there exists a sizeable 
body of literature on phenomenological constitutive equations for the rate- and temperature- 
dependent plastic deformation of metallic materials, to date rational and thermodynamically 
consistent elastic-thermoviscoplastic constitutive relations capable of incorporating the effects of 
large strains and rotations have not been demonstrated. 

Constitutive models for small strain in engineering literature may generally be grouped into 
three categories: classical plasticity, nonlinear visoelasticity, and theories based on microstructural 
phenomena. Each group can be further separated into "unified" and "uncoupled theories, where 
the two differ in their approach to the treatment of rate-independent and rate-dependent inelastic 
deformation. The uncoupled theories decompose the inelastic strain rate into a time-independent 
plastic strain rate and a time-dependent creep rate with independent constitutive relations describing 
plastic and creep behavior. Such uncoupling of the strain components provides for simpler 
theories to be developed but precludes any creep-plasticity interaction. Recognizing that cyclic 
plasticity, creep and recovery are not independent phenomena but rather are very interdependent, a 
number of "unified" models for inherently time-dependent nonelastic deformation have been 
developed recently. 
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Classical incremental plasticity theories are macrophenomenological because they base the 
derivation of state variables purely on experimental results without direct reference to the 
microstructure of the material. Most incremental plasticity theories have four major components: 
(1) stress-elastic strain relations, (2) a yield function describing the onset of plastic deformation, 
(3) a hardening rule which prescribes the strain-hardening of the material and the modification of 
the yield surface during plastic flow, and (4) a flow rule which defines the components of strain 
that are plastic or nonrecoverable. 

Research in this area is voluminous. For example, Zienkiewicz and Cormeau [1] 

< 1 

developed a rate-dependent unified theory which allows for nonassociative plasticity and strain 
softening, but does not model the Bauschinger effect or temperature dependence. Extensions of 
classical plasticity to model both rate and temperature effects were presented recently by Allen and 
Haisler [2], Haisler and Cronenworth [3], and Yamada and Sakurai [4]. 

In the nonlinear viscoelastic approach, the constitutive relation is expressed as a single 
integral or convoluted form. This type of constitutive model employs the thermodynamic laws 
along with physical constraints to complete the formulation. A detailed review of several existing 
theories is presented by Walker [5]. Walker's [5] theory is based on a unified viscoplastic integral 
developed by modifying the constitutive relations for a linear three parameter viscoelastic solid. 
The theory contains clearly defined material parameters, a rate-dependent equilibrium stress, and a 
proposed multiaxial model. An important shortcoming of Walker's theory is its failure to model 
transient temperature conditions. Many other nonlinear viscoelastic theories have been proposed 
including those by Cernocky and Krempl [6], Valanis [7], and Chabache[8]. 

The microphenomenological theories attempt to represent the response of polycrystalline 
materials in terms of various micromechanisms of deformation and failure. Various dislocation 
theories have been developed to predict plastic deformation in terms of dislocation interaction, slip, 
glide, density, etc. Most of the material models developed, to date, depend primarily on the 
number of state variables used and their growth or evolutionary laws. Many of the recent 'unified' 
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microphenomenological theories have been discussed and evaluated by Walker [9] and Chan et al. 

[ 10 ]. 

One example of a microphysically based constitutive law is an elastic-viscoplastic theory 
based on two internal state variables as proposed by Bodner, et al. [11]. These authors, 
demonstrate the ability of the constitutive equations to represent the principal features of cyclic 
loading behavior including softening upon stress reversal, cyclic hardening or softening, cyclic 
saturation, cyclic relaxation, and cyclic creep. One limitation of the formulation though is that the 
computed stress-strain curves are independent of the strain amplitude and therefore too "flat" or 
"square". 

Miller [12] has reported research on the modeling of cyclic plasticity with "unified" 
constitutive equations. He also recognizes the shortcomings of many theories in predicting 
hysteresis loops, which are oversquare in comparison to observed experimental behavior. 
Improvement is accomplished by making the kinematic work-hardening coefficient depend on the 
back stress and the sign of the nonelastic strain term. Theories that are similar in format to Miller's 
have been proposed by Krieg, Swearengen and Rhode [13] and by Hart [14]. The models use two 
internal state variables to reflect current microstructure state and are based upon models for 
dislocation processes in pure metals. All these constitutive theories were formulated without the 
use of a yield criterion. Since these models do not contain a completely elastic regime, the function 
that describes the inelastic strain rate should be such that the inelastic strain rate is very small for 
low stress levels. Theories with a yield function and a full elastic regime have been developed for 
the case of isotropic and directional hardening by Lee and Zavrel [15]. 

As previously noted, the quantities utilized in the small strain theory of viscoplasticity 
(stress, strain, stress rate, and strain rate) are defined only within the assumption of "small strain", 
but this is always left unstated. Whether or not the strains for a given case are "small" cannot be 
determined a priori by geometric considerations. In general, one cannot know in advance whether 
for a given loading of a material the "small strain" assumption (always left undefined) will hold or 
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not. The question of whether the small-strain approximations are valid is always avoided in the 
"small strain" literature. Furthermore, as Hill [16] points out, the really typical plastic problems 
involve changes in geometry that cannot be disregarded. In many cases, for example, it is 
sufficient to take into account finite plastic strains and small elastic strains or vice versa. From the 
theoretical viewpoint it is desirable in all cases to have a theory which intrinsically allows for both 
the elastic and plastic strains to be large. Such a theory of course, must reduce to the earlier 
mentioned special cases, as limiting cases. Furthermore, such theories provide a check for those 
which are obtained by generalizing small strain theories. 

The mathematical theories of deformation and flow of matter deal essentially with the gross 
properties of a medium. Heat and mechanical work are considered as additional causes for a 
change of the state of the medium. The resulting phenomena in any particular material are not 
unrelated. Therefore, a thermodynamical treatment of the foundation of the theory of flow and 
deformation is appropriate, and indeed the obvious approach. Two very different main approaches 
to a thermodynamic theory of a continuum can be identified. These differ from each other in the 
fundamental postulates upon which the theories are based. An essential controversy can be traced 
through the whole discussion of the thermodynamic aspects of continuum mechanics. None of 
these approaches is concerned with the atomic structure of the material. They, therefore, represent 
purely phenomenological approximations. Both theories are characterized by the same 
fundamental requirement that the results should be obtained without having recourse to statical or 
kinetic methods. 

Within each of these approaches there are two distinct methods of describing history and 
dissipative effects: the functional theory [18], in which all dependent variables are assumed to 
depend on the entire history of the independent variables, and the internal variable approach [19], 
wherein history dependence is postulated to appear implicitly in a set of internal variables. For 
experimental as well as analytical reasons[20, 21] the use of internal variables in modeling inelastic 
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solids is gaining widespread usage, in current research. The main differences among the various 
modem theories lie in the choice of these internal variables. 

The predictive value of an elastic- viscoplastic material model for non-isothermal, large 
deformation analyses depends therefore on three basic elements: 

a) the nonlinear kinematic description of the elastic-plastic deformation. 

b) termodynamic considerations, and 

c) the choice of external and internal thermodynamic variables as well as on their interactions. 

The problem of viscoplastic deformations in shells has been treated at several levels of 
approximation and generality. 

The simplest approaches (see [22]) are based on the assumption of infinitesimal 
displacement gradients (which implies infinitesmal strains) and a material model of stationary 
creep, sometimes with an approximate inclusion of primary creep. 

A more general analysis utilizes shell kinematics for moderately large displacement 
gradients (at least some of them), infinitesimal strains, and material models of stationary or simple 
non-stationary creep (see [22]). This type of assumption is capable of solving problems of creep 
buckling [23], and it does reproduce the sometimes stiffening effect of the interaction between the 
normal forces and the normal deflection. Extension of this kind of formulation with a viscoplastic 
material model is presented in [24-26]. The use of numerical methods [27] makes possible the 
solution for many non-trivial types of structures. 

The problems of large strains, which arise in the analysis of large creep or thermal 
deformation of shells, have not been treated at all in a general manner. Recognizing that finite 
strain effects are present in these problems, reliable prediction demands that such effects be 
included rationally and properly in the analysis. In addition to the necessary basic kinematical and 
dynamical equations of the shell theory, such an analysis must incorporate a correctly invariant 
formulation of the material equations and requires an evaluation of the strain-rate tensors through 
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the thickness of the shell. Such an analysis cannot be found in explicit form, at least in the readily 
accessible engineering literature. 

Several authors have developed mathematical description of the kinematics of the three 
dimensional deformation of elastic or viscoelastoplastic materials [28, 29]. However, it is not 
clear how to best select the reference space and configuration for the stress tensor, bearing in mind 
the rheologies of realistic materials. Although an intrinsic relation, which satisfies material 
objectivity can be used [30, 31], the choice is not unique (see for example [29], [32], [33]) . 

For shell-like structures, failure may be caused by buckling. Thus, stability can be a 
primary consideration in the structure design. In a high temperature environment, there will be 
much more concern on this issue, because the inelastic deformation may lead to a geometrical 
imperfection which, in turn, may further decrease the load carrying capability. Therefore, stability 
of shell-like structures can become the main concern for designers, and buckling and postbuckling 
behavior of shell- like structures (see [34]) must be studied. 

In the analysis of shell-like structures, it is worthwhile to note that Donnell [35] and 
Sanders [36] made great contributions in the formulation of nonlinear shell theories. Many 
applications are based on their formulation and simplifications. 

For a long time, the research efforts have been put into the subject on how to improve the 
discrepancy between the theoretical and experimental results. The present trend in buckling and 
postbuckling analyses is to relax several of the assumptions in classical theories and employ 
nonlinear kinematic and constitutive relations. 

As the complexity of shell-like structures increases and as the computational capability 
improves, numerical methods play a major role in predicting buckling and in enhancing our 
understanding of postbuckling response. 

The finite element analysis for shell structures has been the subject of interest for many 
years. There have been many published works in this field. Most of them deal with elastic-plastic 
material behavior. Some of them also deal with geometric nonlinearity and postbuckling behavior. 
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Ahmad, Irons and Zienkiewicz [37] first proposed the 'degenerated' three dimensional 
isotropic element which can be used in the linear analysis of thin shells. The basic assumption 
used in [37] is that the normal will remain normal to the deflected midsurface, straight and 
inextensional after deformation. With this assumption, the displacement field in the shell can be 
expressed in terms of five degrees of freedom (three translations and two rotations) of the nodes 
which are located in the middle surface. 

The 'degenerated' 3-D element was extended to the geometrically nonlinear analysis of 
shells by Ramm [38]. He used both quadratic and cubic interpolation functions in his work. The 

development is based on total Lagrangian formulation. 

Bolouchi [39] developed various shell elements with 8-16 nodes. His work is based on 
both total Lagrangian and updated Lagrangian formulations. 

At the same time, finite element analysis has been adopted to the area of nonlinear 
continuum mechanics. It made it possible to obtain solutions to a large class of nonlinear problems 
with acceptable accuracy. 

J.T. Oden [40, 41] extended the elastic theory of finite elements to the hyper-elastic and 
visco-elastic field. 

In 1968, H.D. Hibbitt, D.V. Marcal and J.R. Rice [42] first introduced hypo-elastic 
formulation into the finite element analysis. They adopted the incremental theory based on a 
Lagrangian reference frame and their formulation is suitable for large strain and large displacement 
response. In their work, they used a linear relation between the Jaumann stress increments and 
increments of the deformation tensor which are invariant with respect to rigid body rotation. The 
formulation for small strain and large rotation approximation is also proposed in their work. 

Needleman [43] derived similar finite element equations from variational principles given 
by Hill [44]. The work is also based on the Lagrangian formulation. 

The Eulerian formulation which is based on the current configuration has been used by 
Yaghmai and Popov [45]. They derived the equations by means of variational principles and 
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solved some elastic and elastic-plastic problems. Similar work was done by Gunasekera and 
Alexander [46]. They gave their equations for Prandtl-Reuss plastic materials and derived the 
equations from the principle of virtual work. 

In 1974, R.M. McMeeking and j.R. Rice [47] introduced Eulerian finite element 
formulation for large elastic-plastic flow, which is parallel to the treatment of [42]. The method is 
based on Hill's variational principle for incremental deformations. It is ideally suited for 
isotropically hardening Prandtl-Reuss materials. 

T.Y. Chang and K. Sawamiphakdi [38] adopted the hypoelastic theory to the degenerated 3- 
D isoparametric element for laminated anisotropic shells. The nonlinear geometric stiffness matrix 
which is based on Lagrangian description was developed and some example problems were 
presented. 

Generally, the Newton-Raphson integration method is used to achieve the correct 
equilibrium positions in the nonlinear finite element computations. However, it is no longer 
appropriate to use it in establishing post buckling response. The reason is that the stiffness matrix 
tends to be singular resulting in an increasing number of iterations. Finally, the result will diverge 
at the critical point. In order to overcome these problems and to trace the response beyond the 
critical point, several efficient methods have been developed [48, 49]. 

Bergan [50] introduced the 'current' stiffness parameter’ to guide the equilibrium 
iterations. When the parameter reaches the prescribed value, the execution of iterations stops. At 
that time, the displacement continues to increase until a new parameter value is reached. This 
means that the iterations are suppressed near the critical point. The algorithm of this method is 
simple and the program is easy to develop from the Newton-Raphson method. The drawback is 
that the load step in the neighborhood of the critical point need still be small. 

Argyris [52] and several other researchers introduced the displacement control method. In 
this scheme, a single displacement component is selected as a control parameter and the 
corresponding load level is considered to be unknown. In the initial paper, the symmetry of the 
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stiffness matrix is damaged. The method, then, was improved as a two step method [52, 53], in 
order to preserve the symmetry of the stiffness matrix. The displacement control method is 
relatively widely used since it is both efficient and easy to control. The limitation of this method is 
that it will fail whenever the structure snaps back from one load level to a lower one. 

Riks [54, 55], and Wempner [56] independently introduced the constant-arc-length 
method. The load increment in this method is confined by the equation of arc length. The 
increments of load and displacement vector are mixed. The final equilibrium position is located by 
continuing drawing a normal plane to the new tangent of the load-displacement curve. Generally 
speaking, this method is efficient in the entire load range and can be applied to all possible 
nonlinear structural responses. 

Crisfield [57, 58] further modified Riks' method. He introduced line searches into the 
constant-arc-length algorithm. In addition, he developed a scheme in which a single parameter is 
employed to accelerate the speed of convergence using the line search concept. His work makes 
Riks’ method more efficient and also much easier for programming. 

2. SUMMARY OF WORK 

The progress made and the work performed have been elaborated upon in an interim 
scientific report submitted to the sponsor in late 1986, and in a series of semiannual progress 
reports. The most recent of these is dated October 1989. 

2.1 Traditional Approach 

Following a traditional approach, a method was developed for bounding the response of 
problems of viscoelastic material behavior, based on nonlinear kinematic behavior. Upper and 
lower bounds are established through bounding of the convolution integral of the governing 
nonlinear Vol terra- type integral equation. Details of the method can be found in the first paper of 
Appendix A. 
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Next, a differential (as opposed to integral) methodology for solution of one-dimensional, 
kinematically nonlinear, visoelastic problems was developed. Using the example of an 
eccentrically loaded cantilever beam-column, the results from the differential formulation were 
compared to the results obtained from the integral solution technique. The details of this are 
found in the second paper of Appendix A. This paper also includes a discussion of the various 
factors affecting the numerical accuracy and rate of convergence of the two procedures. 
Moreover, the influences of some "higher order" effects, such as straining along the centroidal 
axis, are also discussed. 

Finally, an analytic study of beams and arches subjected to significant thermal cycling from 
ambient temperatures up to 800°C is presented. In this study, Walker's [9] unified nonlinear 
hereditary type of viscoelastoplastic constitutive law is employed to characterize the time-and 
temperature-dependent properties of a typical aerospace alloy, Hastelloy X. 

The details are given in the third paper of Appendix A. A shorter version of this paper was 
published in the ASME Journal of Engineering Materials and Technology (January 1990 issue). 
The PVP-Vol. 153 version is made part of this report, because it contains more detail. 

2.2 Novel Approach 

A complete true ab-initio rate theory of kinematics and kinetics for continuum and curved 
thin structures, without any restriction on the magnitude of the strains or the deformation, was 
formulated. The time dependence and large strain behavior are incorporated through the 
introduction of the time rate of the metric and curvature in two coordinate System; a fixed (spatial) 
one and a convected (material) coordinate system. The relations between the time derivative and 
the covariant derivatives (gradients) have been developed for curved space and motion, so that the 
velocity components supply the connection between the equations of motion and the time rate of 
change of the metric and curvature tensors. 

The metric tensor (time rate of change) in the convected material coordinate system is 
linearly decomposed into elastic and plastic parts. In this formulation, a yield function is assumed, 
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which is dependent on the rate of change of stress, metric, temperature, and a set of internal 
variables. Moreover, a hypoelastic law was chosen to describe the thermoelastic part of the 
deformation. 

A time and temperature dependent viscoplastic model was formulated in this convected 
material system to account for finite strains and rotations. The history and temperature dependence 
were incorporated through the introduction of internal variables. The choice of these variables, as 
well as their evolution, was motivated by phenomenological thermodynamic considerations. 

The nonisothermal elastic-viscoplastic deformation process was described completely by 
"thermodynamic state" equations. Most investigators (in the area of viscoplasticity) employ plastic 
strains as state variables. Our study shows that, in general, use of plastic strains as state variables 
may lead to inconsistencies with regard to thermodynamic considerations. Furthermore, the 
approach and formulation employed by all previous investigators lead to the condition that all 
plastic work is completely dissipated. This, however, is in contradiction with experimental 
evidence, from which it emerges that part of the plastic work is used for producing residual 
stresses in the lattice, which, when phenomenologically considered, causes hardening. Both 
limitations are not present in our deformation, because of the inclusion of the "thermodynamic 
state" equations. 

The obtained complete rate field equations consist of the principles of the rate of the virtual 
power and the rate of conservation of energy, of the constitutive relations, and of boundary and 
initial conditions. These formulations provide a sound basis for the formulation of the adopted 
finite element solution procedures. 

The derived shell theory, in the least restricted form, before any simplifying assumptions 
are imposed, has the following desirable features: 

(a) The two-dimensional, impulse-integral form of the equations of motion and the Second 

Law of Thermodynamics (Clausius-Duhem inequality) for a shell follow naturally and 

exactly from their three-dimensional counterparts. 
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(b) Unique and concrete definitions of shell variables such as stress resultants and couples, rate 
of deformation, spin and entropy resultants can be obtained in terms of weighted integrals 
of the three-dimensional quantities through the thickness. 

(c) There are dq series expansions in the thickness direction. 

(d) There is no need for making use of the Kirchhoff Hypotheses in the kinematics. 

(e) Ali approximations can be postponed until the descretization process of the integral forms 
of the First Law of Thermodynamics. 

(f) A by-product of the descent from three-dimensional theory is that the two-dimensional 
temperature field (that emerges) is not a through-the-thickness average, but a true point by 
point distribution. This is contrary to what one finds in the literature concerning thermal 
stresses in the shell. 

To develop geometrically nonlinear, doubly curved finite shell elements the basic equations 
of nonlinear shell theories have to be transferred into the finite element model. As these equations 
in general are written in tensor notation, their implementation into the finite element matrix 
formulation requires considerable effort. 

The nonlinear element matrices are directly derived from the incrementally formulated 
nonlinear shell equations, by using a tensor-oriented procedure. The classical thin shell theory 
based on the Kirchoff-Love hypotheses (Formulation D) was employed for this purpose. For this 
formulation, we are using the "natural" degrees of freedom per midsurface shell node: three 
incremental velocities and the rates of rotation about the material coordinates in a mixed form. 

The quasi-linear nature of the principle of the rate of virtual power suggests the adoption of 
an incremental approach to numerical integration with respect to time. The availability of the field 
formulation provides assurance of the completeness of the incremental equations and allows the 
use of any convenient procedure for spatial integration over the domain V. In the present instance, 
the choice has been made in favor of a simple first order expansion in time for the construction of 
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incremental solutions from the results of finite element spatial integration of the governing 
equations. 

The procedure employed permits the rates of the field formulation to be interpreted as 
increments in the numerical solution. This is particularly convenient for the construction of 
incremental boundary condition histories. 

Even under the condition of static external loads and slowly growing creep effects, the 
presence of snap-through buckling makes the inertia effects significant. In dynamic analyses, the 
applied body forces include inertia forces. Assuming that the mass of the body considered is 
preserved, the mass matrix can be evaluated prior to the time integration using the initial 
configuration. 

Finite element solution of any boundary-value problem involves the solution of the 
equilibrium equations (global) together with the constitutive equations (local). Both sets of 
equations are solved simultaneously in a step by step manner. The incremental form of the global 
and local equations can be achieved by taking the integration over the incremental time step 
t=tj+i-tj. The rectangular rule has been applied to execute the resulting time integration. 

Clearly, the numerical solution involves iteration. A simplified version of the Riks- 
Wempner constant- arch-length method has been utilized. This iteration procedure which is a 
generalization of the displacement control method also allows to trace the nonlinear response 
beyond bifurcation points. In contrast to the conventional Newton-Raphson techniques, the 
iteration of the method takes place in the velocity and load rate space. The load step of the first 
solution in each increment is limited by controlling the length ds of the tangent. Either the length is 
kept constant in each step or it is adapted to the characteristics of the solution. In each step the 
triangular-size stiffness matrix has to be checked for negative diagonal terms, indicating that a 
critical point is reached. 

The study is limited only to the simplest of the developed shell theory formulations 
(Formulation D). 
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The various shell theory approximations (formulations) are based on the use of certain 
simplifying assumptions regarding die geometry and kinematics of the shell configuration. 

These are: 

Assumption I : The material points which are on the normal to the reference surface before 
deformation will be on the same normal after deformation. 

Assumption II : The shell is sufficiently thin so that we can assume linear dependence of the 
position of any material point (in the deformed state) to the normal (to the reference surface) 
coordinate (in the deformed state). The linear dependence can easily be changed to parabolic, 
cubic, or any desired degree of approximation. 

Assumption III : The rate of change of the velocity gradients with respect to in-plan coordinates on 
the two bounding shell surfaces is negligibly small. 

Assumption IV : The rate of change of the distance of a material from the reference surface is 
negligibly small. 

On the basis of the above four simplifying assumptions, several formulations result, for the 
analysis of thin shells. These formulations are denoted below by capital letters. 

Formulation A : This formulation makes use of Assumption I, only. 

Formulation B : This formulation employs Assumptions I and II. 

Formulation C : This formulation employs Assumption I, II and HI. 

Formulation D : This is the classical thin shell theory based on the Kirchoff-Love 
hypotheses of Assumptions I, II, HI, IV, as applied to large deformation theory. 

These formulations are arranged in such a manner that we move from the least restrictive 
(A) to the most restrictive (D). 

In addition to this, a fifth formulation (E) can easily be devised and this formulation in 
terms of order of restriction is similar to Formulation A. Formulation E makes use of Assumption 
II only. 
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Three papers are included as Appendix B. These papers reflect the work associated with 
the novel approach and provide detail in formulation as well as in application. 

3. PUBLICATIONS AND PRESENTATIONS 
Several presentations and publications resulted from this project. Moreover, two Ph.D. 
students were supported (one has already received the degree, the second is expected to complete 
the requirements by the end of 1990), as well as one postdoctoral fellow and three faculty members 
(Drs. R.L. Carlson, R. Riff and G J. Simitses) of the Georgia Institute of Technology. 

A list of the presentations and publications is given below: 

3.1 Presentations 

1. "Thermodynamically Consistent. Constitutive Equations for Nonisothermal, Large Strain, 
Elasto-Plastic Creep Behavior," 26th AIAA/ASME/ASCE/AHS Structures Structural 
Dynamics and Materials Conference, Orlando, FL., April 14-17, 1985. 

2. "Dynamic Creep Buckling: Analysis of Shell Structures Subjected to Time-Dependent 
Mechanical and Thermal Loading," NASA Conference on Structural Integrity and Durability 
of Reusable Space Propulsion Systems, Cleveland, OH., June 4-5, 1985. 

3. "Bounding Solutions of Geometrically Nonlinear Viscoelastic Problems," 27th Structures, 
Structural Dynamics and Materials Conference, San Antonio, TX, May 4-6, 1986. 

4. "Dynamic Analysis of Shell Structures Subjected to Mechanical and Thermal Loading," 
AFOSR/ARO Conference on Non-Linear Vibrations, Stability, and Dynamics of Structures 
and Mechanisms, Blacksburg, VA, March 23-25, 1987. 

5. "Solution Methods for One-Dimensional Nonlinear Viscoelastic Problems," 28th 
AIAA/ASME/ASCE/AHS Structures, Structural Dynamics and Materials Conference, 
Monterey, CA, April 6-8, 1987. 
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6. "Non-Isothermal Elastoviscoplastic Snap-Through and Creep Buckling of Shallow Arches," 
28th AIAA/ASME/ASCB/AHS Structures, Structural Dynamics and Materials Conference, 
Monterey, CA, April 6-8, 1987. 

7. "Creep Analysis of Beams and Arches Based on a Hereditary Visco-Elastic-Plastic Law," 
ASME Winter Annual Meeting, Chicago, IL., Nov. 28-Dec. 2, 1988. 

8. "Non-Isothermal Buckling Behavior of Viscoplastic Shell Structures," 30th 
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Bounding Solutions of Geometrically 
Nonlinear Viscoelastic Problems 

John M. Stubstad* and George J. Similsest 

Georgia Institute of Technology , Atlanta , Georgia 

A method is presented for bounding solutions to problems of linear viscoelastic material behavior formulated 
using nonlinear kinematic measures of deformation. Upper and lower bounds are established through bounding 
of the convolution integral of the governing nonlinear Volterra-type integral equation. A significant feature of 
the manner in which these bounding solutions are generated is that time may be treated as a parameter, thereby 
casting the bounding solutions into a quasielastic context. Consequently, numerical evaluation is simplified since 
the necessity of continually approximating convolution integrals of the deformation history, required for exact 
solution, is eliminated. This, in turn, results in a substantial reduction in the computational effort required for 
numerical evaluation. In one of the example problems considered, this reduction translates into more than a 
thirtyfold difference in computer tiqie needed for determination of the exact and bounding solutions. Applica- 
tion of the bounding technique is demonstrated through two examples and includes a limited comparison with 
some recently published experimental data. 


Nomenclature 

a = load eccentricity 

E\,E 2 = elastic constants of ideal viscoelastic material 

g(s,u) = Green's function for the spatial integrals 

/ = moment of inertia 

J(t) = creep compliance 

$(p) = Laplace transform of creep compliance 

k(t — r) = generic kernel of convolution integral 

K(t) -integrated form of k{t-r) 

L = length of beam 

£ p = Laplace transform operator 

M(s\t) = bending moment at position s' and time t 
p = Laplace transform variable 

P = time-independent load 

P E - Euler load 

R(t) -time-dependent load 

s',s = dimensional and nondimensional distance 

along the beam, respectively 
t,T =time 

x,y - spatial coordinates 

a(r) = angle of rotation of the end of the cantilever 

column 

A(r),6(r) = lateral and transverse deflection of the end of 

the beam, respectively 

Vi = viscous coefficient of ideal viscoelastic material 

= generic representation of a spatial integral 
K(s',t) ^curvature at position s' and time t 

X - scalar parameter 

M = relaxation constant of ideal viscoelastic 

material 

<p(s',t) - angle of rotation of the cantilever at position s' 

and time r 

$(s',p) = Laplace transform of angle of rotation 
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Subscripts 

qe = quasielastic solution 

ubjb = upper and lower bounds, respectively 


Introduction 

I NTEGRAL transform techniques, such as the Laplace 
transform, provide simple and direct methods for solving 
viscoelastic problems formulated within a context of linear 
material response and using linear measures for deformation. 
Application of the transform operator reduces the governing 
linear integrodifferential equations to a set of algebraic rela- 
tions between the transforms of the unknown functions, the 
viscoelastic operators, and the initial and boundary condi- 
tions. Inversion, either directly or through the use of the ap- 
propriate convolution theorem, provides the time domain 
response, once the unknown functions have been expressed in 
terms of sums, products, or ratios of known transforms. 
When exact inversion is not possible, approximate techniques, 
such as suggested by Schapery, 1 may provide accurate results. 

The overall problem becomes substantially more complex 
when nonlinear effects must be included. We consider here 
situations where a linear material constitutive law can still be 
productively employed, but where the magnitude of the 
resulting time-dependent deformations warrants the use of a 
nonlinear kinematic analysis. The governing equations will 
be nonlinear integrodifferential equations for this class of 
problems. Thus, traditional as well as approximate tech- 
niques, such as cited above, cannot be employed since the 
transform of a nonlinear function is not explicitly expressible. 

Rogers and Lee 2 considered such a problem in an in- 
vestigation of the finite deflection of a viscoelastic cantilever 
beam. Employing an analogy of an associated elastic prob- 
lem, they derived a solution to the viscoelastic problem in a 
form involving a time convolution of a nonlinear space and 
time-dependent integral function. Numerical evaluation was 
accomplished using Picard’s method of successive sub- 
stitutions. New'ton-Coates quadratures were employed to ap- 
proximate the spatially dependent integral relationship; a 
mean-value-based finite-difference formula was used for the 
time convolution. 

Solution procedures of this type are generally well suited 
for computer implementation. However, they can become 
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ejcungutataianall^ inefficient w hen? the nespansR must* he deter- 
mined over an extended time period. Each increment in time 
requires a reevaluation of the convolution integrals. Thus, 
the entire deformation history must be retained in memory 
during the calculations. Since each completed set of com- 
putations adds another set of results to this history, this 
generates an ever-increasing memory requirment. In addi- 
tion, the total number of computations performed during 
each succeeding iteration also increases. 

In this regard, an approximation technique proposed by 
Schapery 3 can provide an attractive alternative. Commonly 
referred to as the quasielastic approximation, it has most 
recently been employed by Vinogradov 4 and Vinogradov and 
Wijeweera 5 in studies of the behavior of eccentrically loaded 
viscoelastic cantilever columns. 

The method is based on the observation that the solution 
procedure developed by Rogers and Lee 2 may be interpreted 
as a sequence of short-time interval quasielastic solutions. 
This suggests that approximate solutions may be generated 
by replacing the original viscoelastic problem by an 
“equivalent” time-dependent elastic one. In this replacement 
problem, the elastic properties are equated to the instan- 
taneous values of the relaxation moduli or creep compliances 
of the viscoelastic material. 

The inherent numerical advantage provided by this tech- 
nique is that it eliminates the potentially inefficient calcula- 
tion of convolution integrals. Thus, the speed and efficiency 
at which the time-dependent response is calculated is in- 
dependent of elapsed time. The obvious potential disadvan- 
tage is that, since it is an approximation, significant dif- 
ferences may exist between the actual response of the visco- 
elastic body and those predicted quasielastically. In addition, 
the quasielastic method does not provide a direct method for 
assessing whether any errors incurred are conservative. 

In this paper, we present an approximation technique that 
provides both upper- and lower-bound predictions for the 
class of viscoelastic problems under consideration. From 
these bounds, one may readily deduce when the approxima- 
tion provides sufficiently accurate results or when more in- 
volved techniques must be used. Finally, we demonstrate 
that solutions for this class of viscoelastic problems can be 
accomplished within a Laplace transform context, even 
though the transformed functions cannot be expresssed as 
explicit functions of the transform variable. 

Preliminaries 

As a motivation for the development, consider an integral 
equation of the form 

<p(x>t) = X j o k(t-T)d(X,T)&T (1) 

where <p(x,r) and 0(x,r) may be scalar, vector, or tensor 
functions of position vector x and time r. We shall assume 
that the kernel k(r) is positive semidefinite over the range of 
integration and X is some scalar parameter. In addition, we 
assume that k(r ), <p(x,t), 0(x,r), and their first partial 
derivatives with respect to r are continuous over the interval 
0 + :£r<f. Finally, we assume that y>(x,r) and 0(x,r) are 
continuous with respect to x over some closed domain D 
and possess continuous first and second partial derivatives 
with respect to x over at least the interior of the domain. 

For the class of problems under consideration, the func- 
tion 0(x,t) represents a spatial integral in which <p{x t r) ap- 
pears in the integrand in a nonlinear manner. Depending 
upon the boundary conditions, 0(x,r) may also include ad- 
ditive nonlinear forms of <p(x,t). Thus, Eq. (1) may be viewed 
as a Volterra-type integral equation of the second kind. 

Let us assume that, even before specific solutions have 
been generated, we are able to infer some information about 
nr^neT’Tl manner in which O(x.t) behaves. Suppose, for 


example, that! represents sourer measure; of? dfefIemio.rn 

which (we deduce) must be a nondecreasihg function with 
respect to time. Thus, over the interval 0+ this would 

imply 

0(* > O^<0(*,r)<;0(.M) (2) 

This suggests that if we establish the approximate solutions 


^ = X] o Ar(f-T)0(*.O + )dr 

(3a) 

<fi ub =\^k(t-T) 6 (x,t)dT 

(3b) 


where subscripts ub and lb denote upper and lower bounds, 
respectively, we can then define difference functions A<p ]b 


and A <p ub by 

&<Pib = <p( x pt)~ <Pib ( 4a ) 

&Vub^Vub~<P( X >D ( 4b > 

Then substitution of Eqs. (1) and (3) into Eq. (4) yields 

A^ = \jV(/-r)[0(JC.T)-0(jf,0 + )]dr (5a) 

A?,,, = x( o *(/-r)(0(.*,O-0(*»T)]dr ( 5b ) 


As a direct consequence of Eq. (2), the quantities enclosed 
by square brackets in Eqs. (5) must be positive semidefinite 
for all values of time r. Because both 0(x,O + ) and 0(x,/) are 
constant with respect to r, the square bracket terms must be 
continuous in r since, by prior assumption, 0(x,r) is con- 
tinuous in r. Thus, for a continuous and positive semdefinite 
kernel, the integrand is continuous and positive semidefinite 
over the entire range of integration. Consequently, for 
positive X, the differences A <p lb and A <p ub must be positive 
for all time. Therefore, the approximate solution <p ub must 
represent an upper bound for the exact solution. Similarly, 
(fi !b must inherently bound the exact solution from below. 

The numerical advantages provided by working with the 
bounding functions are easily demonstrated. Letting 

K(t) = xj^ k(t—T)dr (6) 

and noting that 0(x,O + ) and 0(x,t) are independent of r im- 
ply that Eqs. (3) have the form 

* tb =KU)d{x, 0 + ) (7a) 

and 

<p ub =K(t)e(x,t) (7b) 

Thus, the time convolution of the exact solution (Eq. (1)] 
has, in the bounding formulation, been replaced by a format 
in which time appears only as a parameter. Consequently, 
numerical solution of Eqs. (7) requires integration only over 
the spatial coordinates, whereas the exact solution requires 
both spatial and temporal integrations. 

The preceding development was based upon the assump- 
tion that S(x,t) was a nondecreasing function with respect to 
time. The technique is easily adapted to cases where 0(x,r) is 
a nonincreasing function. Thus, if 
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9(x,0*)z9(x,r)^9(x,t) 


( 8 ) 



we can replace Eqs. (3) with 


V>/» = xj o k(t-r)$(x.t)dr 

(9a) 

^ = xj'*(/-r)0(x,O + )dT 

(9b) 


Thereafter, proceeding as before generates the desired 
bounding behavior. In a similar manner, a simple series of 
modifications to the definitions of the bounding functions 
are needed when X is a negative rather than positive scalar. 

Applications of the bounding technique, including com- 
parisons with exact solutions, are provided in the following 
sections. 


Applications 
End-Loaded Cantilever Beam 

As noted earlier, Rogers and Lee 2 developed the first, and 
a numerically exact, solution for the problem of an end- 
loaded linearly viscoelastic cantilever beam. The solution, in 
the general form of Eq. (1), was evaluated by employing 
Newton-Coates quadratures for the spatially dependent in- 
tegral function 6(x,t) and a mean value based finite- 
difference formula for the convolution with k(t — r). Several 
years later, Schapery, 3 in a paper describing the quasielastic 
method, presented an approximate solution for this problem. 
Since this approximate solution was in the form of Eqs. (7), 
Schapery was able to generate numerical results directly from 
the elastic analysis presented in Ref. 2. Thus, his solution re- 
quired only numerical evaluation of a spatially dependent in- 
tegral equation with time treated as a parameter. 

Here, we analyze the same problem using the bounding 
procedure. It is demonstrated that Schapery’s approximate 
solution is, in fact, a lower-bound solution for a suitably 
restricted range of deformation. In addition, it is shown that 
a reasonably close upper-bound solution may be readily ob- 
tained. 

Derivation of the governing integrodifferential equation is 
documented in Ref. 2 and thus only summarized here. The 
beam is assumed to be thin and composed of a linearly visco- 
elastic material. Its geometry in the deformed configuration 
is illustrated in Fig. 1. The loading is assumed to be applied 
quasistatically and thus inertia terms are neglected. 

Reference line extensional strains are assumed to be 
negligibly small. Thus, a coordinate s' is employed to specify 
position in both the initial and the deformed states. A non- 
dimensional coordinate s is defined by dividing s' by the 
beam length L . Assuming a linear distribution of strains 
through the depth, bending thus occuring within a Bernoulli- 
Euler context, results in the moment-curvature relationship 
given by 

(t)L R r% <■») 


dx(s\Q 

ds' 


cos 


dyjs\t) 

ds' 


sin 


(12b) 

(12c) 


Substitution of Eqs. (11) and (12a) into Eq. (10), differentia- 
tion with respect to s’ and use of Eq. (12b) yield, after non- 
dimensionalization 


d 2 *(s,t) 

ds 2 


£)i>-«r < T to > 


The associated boundary conditions are 


(13) 


v><0,/)=0 


(14a) 


Mht) 

ds 


(14b) 


It is assumed the beam is undeformed for r< 0. Then, by 
taking the Laplace transform of Eqs. (13) and (14), we 
obtain 


d 2 $(s t p) 

ds 2 

= - (-L^ /> f|(p)je p (/{( T )cos^(5,r)] 

(15a) 


*( 0,p)=0 

(15b) 


d*(l./>)_ 0 

ds 

05c) 

where £ p [ ] 

is the Laplace transform operator, 

p the 


transform variable, and $(p) and 4>(s,/?) the transforms of 
J(t) and <p(s,t), respectively. Here, we have tacitly assumed 
that the Laplace transform of the expression within the 
brackets exists, in the usual sense, even though a formal ex- 
pression for it is not available. 

Assuming that the transform variable appears only as a 
parameter, Eq. (15a) can be viewed as a type of ordinary dif- 
ferential equation. Integration thus yields 

34>(5,p) a$>(0 ,p) 

ds ds 

- _ ^^jpS(p) j o e-^ T /?(r)cos^(w,r)drjdw (16) 

where the Laplace transform has been expressed in the ex- 
plicit manner. Application of the boundary condition, as 


where *($',/) denotes the curvature and A f(s',r) the bend- 
ing moment at location s'. I is the moment of inertia of the 
beam and J(t) the creep compliance of the material. The 
moment at position s' is given by 

M(s’.r) =R(t) [L -x(s',t) - A(r)] (11) 


where R(t) is the end load (see Fig. 1). 

From kinematic considerations, we note that 




dy(s',t) 

ds' 


(12a) 


L 
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Fig. 1 Geometry of the end-loaded cantilever beam. 




0 10 20 
TIME, HR 

b) At loaded end. 


Fig. 2 Angle of deflection of the end-loaded cantilever beam. 


given by Eq. (15c), yields 

(’t) /,5(p) 1o" e ~‘ pT/? < T> [jj C0S v?(M.r)d*/]dr 

(17) 

Note that we have interchanged the order of integration of 
Eq. (16). This follows directly from the assumption of inex- 
tensionality; therefore, s and r represent independent vari- 
ables. Substitution of Eq. (17) into Eq. (16) yields 


The term in the brackets may be simplified through integra- 
tion by parts and by employing 

g(s,u)~i4 t 0 <W<5 

= 5 , ssu^l ( 20 ) 

to yield 

$(s,p) = (~t)p31p)£ p [rIt) g(5,w)cosy?(w,r)dwJ (21) 
Thus defining 0(s,r) by 

d(5,r) = | o ^(j,w)cosv?(w,r)dw (22) 

results in 

*(s,p)~(L l /I)p$(p)£ p [R{r)e(s t T)] (23) 

Next, the Laplace convolution theorem is employed to invert 
Eq. (23) to yield 

«<->- (t) J.W <*> 

Upon a final integration by parts we have 

>p(s,t) = (^p) [/(O)/?(/)0(s,/) + jV(/-r)/?(r)0(s,r)dr] 

(25) 

where ( )' denotes differentiation with respect to the argu- 
ment of the function. Note that Eq. (25) is the viscoelastic 
solution reported in Ref. 2. 

Bounding solutions are developed in the following man- 
ner. Provided, after quasistatic application, the load R(t) is 
held constant at some value P , it is reasonable to presume 
that the angle of deflection y>(s,r) will be a nondecreasing 
function in time. Thus, for the interval 0 + <rst, 

ip(s,0 + )^v(s,T) (26) 

Consequently, restricting our attention to a range of deflec- 
tion such that 0<v?<7r/2 we may conclude that 

cos*? (s,0 + s cosy? (s, T ) > cosy? (s, t) (27) 

Thus, from Eqs. (22) and (27) we have 

9(s,0 + )l>d(s t T)>:0(s,t) (28) 

Through the use of Eq. (28), we can bound the convolution 
integral of Eq. (25) as follows: 


X cosy?(w,r)duJdr 


P0(s,O + )\ o J'U-T)dTS:^ o J'(t-T)R(T)9(s,T)dT 


Integration of Eq. (18) with respect to s, use of the bound- 
ary condition [Eq. (15b)], and manipulation as before yield 



(f-r)dr 


(29) 


*(s,p) = - 


(~r) p3ip) lo 


e-^R(T) 


Integration of the first and third integrals in Eq. (29) and 
substitution of these results into Eq. (25) yield, after 
rearrangement, 




cosy?(w,r)dudr dr 


— J(t) 
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b) At loaded end. * 

Fig. 3 Vertical deflection of the end-loaded cantilever beam. 

Vub (s,t) =7(0) (~) («<*./) + 1]*(*.0 + )] 

(30b) 

Note that Eq. (30a) is the quasielastic solution proposed by 
Schapery. 

Bounding of the deflection of the beam can be readily ac- 
complished by using Eqs (30). Nondimensionalization and 
then integration of Eq. (12c) yield 

r s in*>(M,OdM (31) 

L Jo 

Thus, since <p lb (s,t)<<p(s,t)<<p ub (s,t), we note that, for 
0< <£>< x/2 

sin tp lb (s,t)< sin</? ( s,t ) £ sin?,* ( s,t ) (32) 

which yield, upon substitution into Eq. (31) 

f’sin^O/.Odu (33a) 

L Jo 

- ub - = ( siny’u;, (M,/)du (33b) 

L Jo 

Numerical solutions of Eqs. (25), (30), and (33) are 
generated in a manner similar to Ref. 2. Picard s method of 
successive substitutions is employed to numerically solve the 
integral equations, Eqs. (25) and (30). Spatial integrals are 
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Fig. 4 Geometry of the eccentrically loaded cantilever column. 



Fig. 5 Ideal three-element “limited*’ creep model. 


approximated by using Newton-Coates formulas; a fixed- 
step trapezoidal rule is used for the time convolution. All 
computations are performed on a CDC Cyber 180/855 
located at the Georgia Institute of Technology. 

Figures 2 and 3 compare the results obtained with the 
bounding formulation to the exact solution for the loading 
case reported in Refs. 2 and 3. The form of the creep com- 
pliance for this example is 

J(t)/J(0) = 1 + 7.6x 10“ 4 / + 1.12(1 -e~ 0055; ) (34) 

Figure 2 demonstrates that the bounding solutions provide 
a reasonably narrow band at both the midspan and the 
loaded end locations. For this particular case, the lower 
bound tends to more closely approximate the actual solution. 
For both the upper and lower bounds, the discrepancy be- 
tween the approximate result and the exact solution tends to 
increase with elapsed time and distance from the clamped 
end. 

Figure 3 compares the calculated vertical deflections at 
midspan and at the loaded end. It can be noted that the 
descrepancies between the bounding and exact solutions 
behave in the same manner as described above. In terms of 
absolute accuracy, after 24 h, the upper-bound end deflec- 
tion exceeds the exact result by approximately 3. 3 Vo. The 
lower bound, at the same point and time, is only 1.2% less 
than the exact deflection. Using a 0.1 h fixed-length time 
step, the exact solution required 15.6s of CPU time to com- 
pute. Calculated for the same number of time intervals, each 
of the bounding solutions required only 4.1 CPU s. Note, 
however, that the accuracy of the bounding solutions is in- 
dependent of the length of the time step. Thus, identical 
bounding solution results can be obtained with, for example, 
a 1.0 h time step. Using this larger time step, the computa- 
tion time for each of the bounding solutions was reduced to 
only 0.5 CPU s. Accurate results could not be obtained from 
the exact solution using a time step this large due to errors in 
approximating the convolution integrals. 

Eccentrically Loaded Cantilever Column 

Although formally similar to the prior example, this prob- 
lem is of interest because the eccentric loading generates ad- 
ditive nonlinear boundary terms in the general solution. The 
presence of these terms substantially influences the accuracy 
of the predictions of the system response. 
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We note that the boundary conditions for this problem are 



a) a/L= 0.01. 



Fig. 6 End deflection of the eccentrically loaded cantilever column 
for the three-element model. 


A quasielastic solution of this problem was recently 
presented by Vinogradov. 4 Included with the analysis were 
numerical results for two ideal constitutive models, using 
two eccentricity ratios and a wide range of applied loads. 
Subsequently, Vinogradov and Wijeweera 5 and Wijeweera 6 
published comparisons of results obtained using the 
quasielastic approximation of Ref. 4 to experimental data 
from tests conducted on PTFE G-700 columns. The loading 
and eccentricity ratios employed in those tests were, 
however, restricted to relatively narrow ranges in value. 

Bounding solutions for the problem are developed in a 
similar manner to the prior example. The column is assumed 
to be inextensional, linearly viscoelastic, and loaded quasi- 
statically. Its geometry in the deformed configuration is il- 
lustrated in Fig. 4. Note that the applied load remains 
parallel to the x axis. 

For this geometry, the moment at any position s' is given 
by 


M(s\t) = /?(r) [ 5 ( r ) -Facosa(r) - y(s',r ) } (35) 


Substitution of Eq. (35) into Eq. (10), differentiation with 
respect to s' , use of Eqs. (12a) and (12c) followed by non- 
dimensionalization yield 


aVteO (L*\ f' y/ , Jd[Rlr)sirup(s t r)] 

S? T, 


ip(0,t) -0 

(37a) 

M(\,t)~aR ( f )c 0 Sy> (1,0 

(37b) 

where, for a rigid “extension,” 


«<0-*(U> 

(37c) 


Through the use of Eqs. (10) and (12a), the second boundary 
condition [Eq. (37b)] can, after nondimensionalization, be 
expressed entirely in terms of p by 


= (.*.) j'__ A ,. r) j-iif*w^U.] df(37 d) 

Assuming that the column is undeformed for r< 0 and 
that following the procedure detailed in Eqs. (15-25) again 
yield a solution of the form of Eq. (25) except that, in this 
case, 

0(s,l) = (y) cos ' p < 1,/ * + (38) 

Note that the nonlinear boundary term cosy>(l,r) appears in- 
side the convolution integral as well as in the integrated term 
of Eq. (25). 

Under a constant load P it is again plausible to assume 
that <p(s,r) will be a nondecreasing function with respect to 
r. Thus, in addition to Eqs. (26) and (27), we note 

sin^ (5,0 + ) < simp (s, r ) < sin^ (s, t ) (39) 

Because of the differences in bounding behavior in Eqs. (27) 
and (39), the convolution integral of the general solution 
[Eq. (25)] is split into two separate integrals that are bounded 
individually. 

Substitution of the appropriate bounding functions from 
Eqs. (27) and (39) into Eq. (38) and substitution of these 
results into Eq. (25) yield, after integration and rearrange- 
ment of terms 


Vubis.t) = •/( 0) (pf) [(x) C0S ^ (l 

7 / *\ 1*1 


i_ i/\V/ ' 


(40a) 


(40b) 


Numerical evaluation of Eqs. (40), as well as the exact solu- 
tion given by Eqs. (25) and (38), is accomplished in the same 
manner as the prior example. Figures 6 and 7 present the 
results of these computations for the ideal “limited creep” 
model used in Refs. 5 and 6 and illustrated on Fig. 5. For 
this particular constitutive model, the creep compliance J(t) 
has the form 






b) a/L = 0.10. 

Fig. 7 End deflection of the eccentrically loaded cantilever column 
for the three-element model. 

where we have employed 

7(0)= \/E 2 (41b) 



13 5 7 


time.hr 

a) P/P E = 0.40 and a/L = 0.07 . 
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b) P/P E = 0.125 and a/L = 0.03. 

Fig. 8 Reported and normalized experimental data for the six 
element model. 


P=E l/l?! ( 41c ) 

Figure 6 presents results for E 2 /E x -(}.5 and <i/L = 0.01 
for a range of load ratios. In this figure, and all succeeding 
ones, we employ a viscoelastic “Euler” load P E to non- 
dimensionalize the loading, where 


x 2 / 

Pe ~ 47(0)I 2 


(42) 


For comparison, this Figure also includes results from a 
“standard” quasielastic solution <p ge of the general form 





(43) 


where the various functions on the right-hand side are as 
previously defined. 

It can be observed that the quasielastic solution is almost 
identical to the upper-bound result for the load ratios of 0.67 
and 0.75. At the load ratio of 0.50, the upper-bound and 
quasielastic results differ only in the fourth decimal 
place. Thus, only the upper-bound result has been indicated 
in the Figure for this load case. 

Although the upper-bound and quasielastic results com- 
pare favorably with each other, neither of them r-r the 


lower-bound solution provides a good approximation to the 
exact result except at the lowest load ratio, 0.50. Thus, 
reliance on only a quasielastic type of solution, especially for 
the higher loading instances, could lead to erroneous results. 

With only a quasielastic solution, it is impossible to deter- 
mine its accuracy without calculating the exact solution. 
Thus, it is not possible to assess the magnitude or character 
(i.e., conservative or nonconservative) of the potential er- 
rors. In contrast, the amount of separation between the 
bounding solutions provides such a capability. The narrow 
separation evident at a load ratio of 0.50 might well provide 
sufficiently accurate results without resorting to the more in- 
volved analysis. The significant differences between the 
bounds at the higher loads, instead, indicate that exact solu- 
tions must be determined for accurate results. 

Figure 7 provides results for the same ratio of modulii, but 
for a load eccentricity of 0.10. Again, at the lowest load 
ratio (0.30), the quasielastic and upper-bound solutions are 
virtually indistinguishable and only the upper bound is 
indicated. 

Comparison of Figs. 6 and 7 illustrates that the increase in 
load eccentricity generates several pronounced effects. The 
quasielastic solution, in general, tends to prpvide a more ac- 
curate prediction of behavior at all load levels for the higher 
load eccentricity. Additionally, the larger eccentricity tends 
to decrease the spread between the upper- and lower-bound 
approximations. This is no:, however, a complex 1- - rv^fr?. 1 
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Fig. 9 End deflection of the eccentrically loaded cantilever column 
for the six-element model. 


trend, since at the load ratio of 0.50 the bounding is much 
more narrow for the lower eccentricity case. 

Since some experimental results also are available, 5,6 a 
comparison with these data is worthwhile. Figure 8 provides 
reported as well as “normalized” data for two load- 
eccentricity cases. The normalized curves are generated by 
adding the relative displacement of a specimen from its 
time = 0 + (i.e., 10 or 20 s deflection) to the average 0 + 
displacement of that test group. In this way, the significant 
differences between the observed results in a test group due 
solely to the differences in “instantaneous” deflection could 
be eliminated. As indicated in the figure, this virtually 
eliminates the substantial differences between observed 
results. 

Based on data from four point bending tests the specimen 
. material (PTFE G-700) was modeled 6 as a six-element 
“unlimited creep” type of material. The numerical form for 
the creep compliance is given by 

J(t)/J(0)= 1 + 3.7x 10~ 4 f + 0.17(1 -<?- 902 ') 

+ 0,13(1 -e” 004 ') (44) 

Figure 9 presents the comparison of the exact quasielastic, 
upper- and lower-bound results to the normalized test data 
of Fig. 8. Note that, while there is an apparent significant 
difference between observed and calculated results for the 
higher-load case, differences in the 0 + deflection account for 
most of it. The average reported “instantaneous” non- 
dimensional deflection was 0.0633, whereas the calculated 


value was only 0.0580. If the various results were to be nor- 
malized to eliminate this difference, the test data band would 
completely overlap the calculated results. However, using 
such a procedure for the lower-load case would decrease the 
correlation indicated on Fig. 8. Since the observed “instan- 
taneous” deflection was only 0.00504, normalizing the data 
to the calculated deflection of 0.00530 would move the band 
of test data so that it would be somewhat above the 
calculated results. The main conclusion to be drawn is that, 
qualitatively, the calculated results agree with the observed 
data. Exact comparability is, however, hindered by the large 
differences in initial displacement evident in the test data. 

Conclusions 

A methodology is presented wherein problems of isother- 
mal linear viscoelastic behavior, formulated using nonlinear 
kinematic measures of deformation, may be analyzed 
through the use of a bounding procedure. The bounding 
solutions developed by this technique are similar in form to 
that of a time-dependent elasticity problem. As such, 
numerical solutions may be generated without requiring the 
computation of convolution integrals of the entire history of 
deformation. In one of the examples considered, it is shown 
that this results in an increase in computational efficiency 
more than 30 times greater by comparison to the more tradi- 
tional approach. 

It is also demonstrated that the bounding procedure pro- 
vides reasonably accurate results for a variety of loading 
conditions. In those cases where narrow bounds cannot be 
established, it is shown that a standard type quasielastic ap- 
proach is not necessarily more reliable. The clear implication 
of the wide bounds is that the more involved traditional ap- 
proach must be employed if highly accurate results are 
required. 

As presented, the bounding technique can be directly 
employed for problems where the governing functions may 
be characterized as either nonincreasing or nondecreasing 
functions with respect to time. However, the range of ap- 
plicability potentially can be expanded to include some forms 
of multimodal functions. In general, this would require that 
these functions be capable of being characterized, at least in 
a piecewise manner, as a sequence of unimodal segments. 
Similar to the procedure that was employed in the second ex- 
ample, each of these segments would then be bounded in- 
dividually. The degree of accuracy that might be obtained 
using such a procedure, however, requires further study. 
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Solution Methods for One-Dimensional Viscoelastic Problems 


John M. Stubstad* and George J. Simitsesf 
Georgia Institute of Technology, Atlanta , Georgia 


A recently developed differential methodology for solution of one-dimensional nonlinear viscoelastic problems is 
presented. Using the example of an eccentrically loaded cantilever beam-column, the results from the differential 
formulation are compared to results obtained from a previously published integral solution technique. It is shown that 
the results from these distinct methodologies exhibit a high degree of correlation with one another. A discussion of 
the various factors affecting the numerical accuracy and rate of convergence of these two procedures is also included. 
Finally, the influences of some “higher-order” effects, such as straining along the centroidal axis, are discussed. 


Nomenclature 

a » load eccentricity 

A = area of cross section 

E = Young’s modulus 

g(j,u) = Green’s function for the spatial integrals 
I — moment of inertia 

j(/) = creep compliance 

/ = length of cantilever 

L p = Laplace transform operator 

M,N — moment and force resultant, respectively 

P - applied load 

P e = Euler load 

s, s = dimensional and nondimensional distance along 

the beam, respectively 

t, x = time 

u, w = axial and transverse displacement, respectively 

x , y = spatial coordinates 

y q = Newton-Cotes quadrature weights 

£ 0 = centroidal axis strain 

k ~~ = curvature 

<f> — angle of rotation 

Introduction 

A NUMBER of solution methods are available for vis- 
coelastic problems in which the behavior of the material 
may be adequately characterized by a linear viscoelastic opera- 
tor and where the deformation of the body is sufficiently small 
to allow the use of a linear kinematic formulation . 1,2 Com- 
monly, integral transform methods, separation of variables, 
series expansions, and other techniques provide methodologies 
wherein exact closed-form solutions may be derived. When 
exact solutions cannot be obtained, approximate techniques, 
such as one proposed by Schapery , 3 provide an alternate 
approach. 

The inclusion of nonlinear effects in the analysis significantly 
reduces the mathematical tractability of the problem. These 
nonlinear influences can be induced by geometric factors re- 
sulting .from the magnitude of the deformation or from gross 
rotation of cross sections. Alternatively, nonlinearities in the 
material response may need to be included to provide an accu- 
rate model for material behavior. 

Independent of whether the nonlinearities are produced by 
geometric or material effects, they invariably result in non- 
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linear governing equations. Thus, the solution methods men- 
tioned here, applicable to linear problems, generally cannot be 
employed. Approximation methods , 4 however, have been de- 
veloped and can be employed to analyze such problems. 

One of these methods is to idealize the problem in a manner 
that inherently simplifies the governing relations. An example 
of this technique was the utilization of an ideal “I" cross- 
sectional geometry in early column creep-buckling studies . 5 
With this approximation, the equilibrium equations were re- 
duced to simpler forms, involving the “average*’ stresses in the 
ideal flanges, where closed-form solution was possible. 

Another approach used extensively was to restrict consider- 
ations to only certain types of time-dependent material behav- 
ior . 6 In some cases, this involved retaining only secondary 
creep behavior in the material model. Alternatively, and espe- 
cially when “power-law” type constitutive laws were used, the 
constants or exponents of the law were restricted to special 
.values for which closed-form solution was possible . 7 In a few 
cases, this approximation, as welLas the aforementioned geo- 
metric simplification technique, were employed simultaneously 
to enable solution. A survey of many of these techniques has 
been provided by Hoff . 8 

A more general technique for the solution of geometrically 
nonlinear viscoelastic problems was first presented by Rogers 
and Lee . 9 In this method, the solution was formulated as an 
integral equation that was nonlinear in both time and space. 
From this, numerical results were obtained using computa- 
tional techniques. A recent paper by the authors 10 provided a 
method for bounding the solution of problems formulated in 
this manner. 

Generally, both the exact and bounding technique can be 
employed for problems wherein the response of the material 
may be adequately characterized using a linear viscoelastic 
model, but where the resulting time-dependent deformation of 
the body warrants the use of a nonlinear kinematic formula- 
tion. Problems involving nonlinear viscoelastic material behav- 
ior, however, currently cannot be addressed with this method. 
Unfortunately, many materials, and especially the elevated 
temperature behavior of most metals, require a nonlinear con- 
stitutive characterization. Consequently, an alternate solution 
procedure for one-dimensional problems involving nonlinear 
kinematic and nonlinear material effects has been developed. 
This method, hereinafter referred to as the differential formula- 
tion, is based on the direct solution of the nonlinear differential 
equations of equilibrium. 

Similar to the integral method, the differentia! formulation is 
predicated on the assumption of a quasistatic response. This, 
effectively, “decouples” the temporal and spatial dependence 
of the problem in a manner that allows the general solution to 
be treated as the sequential combination of solutions to a non- 
linear "boundary value” problem and a nonlinear “initial 
value problem. The first of these, the equations characterizing 
the time-dependent states of quasistatic equilibrium, are solved 
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through the use of a Newton-type method. 11 The "‘initial 
value” problem, resulting from the nonlinear constitutive law. 
governs the manner by which the body progresses from one 
state of quasistatic equilibrium to the succeeding one. Numeri- 
cal solutions for this part of the problem arc generated using a 
fourth-order Rungc-Kutta method. This general method has 
recently been employed to examine the nonlinear thermovis- 
coelastic behavior of thin structural members. 12 

In addition to presenting the differential formulation tech- 
nique, a comparison of results obtained using the integral and 
differential formulations is provided. The problem of an eccen- 
trically loaded viscoelastic cantilever beam-column is employed 
as the vehicle through which the comparison is performed. 
Because of the inherent limitation of the integral technique, 
this comparison is restricted to the consideration of a linear 
viscoelastic material. The specific case considered is that of the 
three-parameter viscoelastic solid, which has been examined in 
a number of studies. 10,13 ,4 The results obtained from these two 
distinct methods of solution exhibit a surprisingly high degree 
of correlation with one another, thereby establishing a high 
level of confidence in the validity of the methods. Finally, the 
differential formulation is employed to examine the influences 
of some “higher-order” effects in the class of problems under 
consideration. 


Integral Formulation 

The Rogers and Lee formulation method, 9 hereafter referred 
to as the integral solution technique, is focused toward formu- 
lating the solution to the nonlinear viscoelastic problem in 
terms of an integral equation. The general method was evolved 
through analogy to the associated geometrically nonlinear elas- 
tic problem. Only a synopsis of the method is presented here 
since complete developments for the technique are available in 
the literature. 9,10 

In the integral formulation, the time dependence of the mate- 
rial behavior is expressed in the form of a Volterra-type inte- 
gral operator. This operator acts upon- a second integral 
expression, which characterizes the quasistatic equilibrium of 
the body. For application to beam-column-type problems, it is 
assumed that the beam-column is thin and composed of a lin- 
early viscoelastic material. In addition, referenced line exten- 
sion strains are assumed to be negligibly small. Thus, the 
coordinate £, denoting distance along the undeformed length, 
can be employed to specify position in both the initial and 
deformed configurations. For convenience, a nondimensional 
coordinate s is defined by dividing s by the length of the beam 
/. Figure 1 illustrates a typical geometry used with this method. 
For the sample problem, the eccentric load is assumed to~be 
applied quasistatically, and its direction does not vary with 
time. 

Assuming a linear distribution of the strains through the 
depth, bending, thus occurring within an Euler-Bernoulli con- 
text, results in a moment-curvature relationship of the form 



J{! - T) |“a^J dt 


(D 


where* k(£,/) denotes the curvature and A/(£,t) the bending mo- 
ment at location s. / denotes the moment of inertia of the beam 
and J{t) the creep compliance of the material. For the sample 
problem, the moment at position s would be given by 


- P(x)[S(x) -h a cosa(t) — y(s, x)] (2) 

where P(x) denotes the load applied at eccentricity a. Since 

d<t>(sj) 




ds 


( 3 ) 


then, for quiescent initial condition, the Laplace transforma- 
tion. of Eq. ( 1) yields 


= -pJ(p)L p [M(s,t)] 
as I H 


( 4 ) 


where }(p) and <1 denote the transforms of J(t) and 
respectively, and where p represents the Laplace variable. Note 
that the nondimensional coordinate s has also been employed 
in the preceding expression. 

Assuming that the Laplace variable apears only alge- 
braically, Eq. (4) has the form of a type of “ordinary” differen- 
tial equation. Consequently, integrating with respect to s yields 


i(s,p) - <i>(0 ,p) = jpJ(p)L p U dr 

0 


( 5 ) 


Note that the order of integration and Laplace transformation 
has been interchanged. This is a direct result of the assumption 
of inextensionality; consequently, s and t represent indepen- 
dent variables. 

Equation (5) reveals a very interesting aspect of this formu- 
lation. Namely, the underlying structure of the equation is 
completely determined by the manner in which the moment 
depends on the deformation. For example, even if the moment 
depends upon the spatial coordinate 5, provided it is indepen- 
dent of the deformation, then the basic equation is, in princi- 
ple, integrable to a closed-form solution. This solution is, in 
fact, the usual result obtained from a linear analysis. 

Illustrating how the equation structure changes when the 
moment is related to the deflection is best demonstrated 
through analogy with the associated elastic problem. Note that 
the governing relation for the associated elastic problem can be 
obtained by replacing the creep compliance by the elastic com- 
pliance and eliminating the Laplace operator. This yields 

WS)-^(0)=^J^A/(r)dr ' ..(6) 

Note that the governing equation for the associated elastic 
problem takes on the form of a linear Fredholm equation when 
the moment depends linearly upon the deflection. In contrast, 
a nonlinear Fredholm format is obtained for cases where the 
moment is nonlinearly dependent upon deformation. In a sim- 
ilar manner, the viscoelastic problem takes on a linear “quasi- 
Fredholm” form when the moment is linearly dependent upon 
the deformation. A nonlinear “quasi-Fredholm” format occurs 
when, as in the sample problem, the relationship between mo- 
ment and deflection is nonlinear. 

To complete the formulation for the sample problem, Eqs. 
(2) and (3) are substituted into Eq. (1). Following differentia- 
tion with respect to j, the kinematic relations 


dx(s y t) 

~sT 


= cos <£($,/) 


Syilt) 

ds 


= sin <£(£,/) 


(7) 
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Fig. I Beam-column geometry for the integral method. 



arc employed to express all deformation-related quantities in 
terms of </>. For the sample problem, the boundary conditions 
are 


</><0./) =0 

M(U) « aP{t) cos</>(/,f) (8) 


where ct(/) = <t>(U) for a “rigid" extension. Thus, using the 
methodology detailed in Ref. 9 and assuming quiescent initial 
conditions yields the solution 


<t>(sj) 



J(0)P(t)O(s,t) + 


I 


JV 


- t)/ > (t)0(5,t) dt J 
(9) 


where 

0 ( 5,0 = cos<p( 1,/) + J g(5,r) sin^(r,/) dr ( 10) 
and 

g(5,r) ~ r, 0^r^5 ^ 

= 5 , s £ r £ 1 

Note that the prime in Eq. (9) denotes differentiation with 
respect to the argument of the function. 

Equation (9) represents a time convolution of a nonlinear 
spatial integral equation, Eq. (10). Numerical solutions are 
obtained using Picard’s method of successive substitutions. 15 
Newton-Coates formulas are used to approximate the spatial 
integral, and a fixed-step trapezoidal rule is employed for the 
time convolution. The general format of the algebraic expres- 
sions obtained in this manner is 

«*i.O = T 7 A {(‘ + \ 

' + "l JVn ~ {,)©(*,.',) + ^ /'(f„)©fe0 + ) j (12) 

with * 

S'O r ' 

0(5,,/,) = -y coscj>(l,tj) + &s\ Y. ViSk sin<f>{r k ,tj) 

1 - . L *-‘ 

+ £ 'M sin<f>(r k ,tj) J (13) 

Note that in Eq. (12), -the number of terms in the summation 
increases linearly with each succeeding time step, whereas the 
number of terms in the 'Summations represented by Eq. (13) is 
fixed. This increasing summation requirement in Eq. (12) has a 
significant impact on the relative speed of the integral formula- 
tion computations. 

Differential Formulation 

As previously noted, the differential formulation technique is 
based bn the direct solution of the governing differential equa- 
tions. Similar to the integral formulation, the differential for- 
mulation is also based on the assumption of quasistatic 
behavior. From this, the equations governing the successive 
states of quasistatic equilibrium may be expressed in terms of 
deformation functions and force and moment resultants. Con- 
sequently, these equations have the gencra]_format of a nonlin- 
ear boundary value problem. A Newton-type method, first 
suggested by Thurston, 11 is employed to derive solutions for 
this part of the problem. 

On the other hand, the constitutive law, expressing the time 
dependence of the material response, governs the evolution of 
the force and moment resultants as the system progresses be- 
tween successive states of quasistatic equilibrium. This repre- 


sents a form of initial value problem with the values of the 
constitutive variables, such as accumulated viscoelastic strain, 
providing the initial conditions. For a nonlinear constitutive 
law, numerical procedures such as a Runge-Kutta or Euler 
method may be employed to predict the growth of these vari- 
ables. Note that, for a “beam-theory" type formulation such as 
this, a spatial integration of the viscoelastic strain across the 
cross section is also required to enable evaluation of the force 
and moment resultants. 

Similar to the integral formulation, the differential formula- 
tion for the sample problem is also based on the assumption 
that bending of the beam occurs in accordance with the Euler- 
Bernoulli hypotheses. Employing the functions u(s 9 t) and \v(sj) 
to denote, respectively, the axial and transverse deflection of 
the centroidal axis, then the extensional strain at the centroidal 
axis e 0 is approximately given by 




(14) 


Note that the term ^ du/ds ) 2 has been neglected as small in 
comparison to cu/ds . If, in addition, both the strain at the 
centroidal axis and cu/vs are small in comparison to 1, then it 
is simple to show that 


d<t> ^ dwd 2 u dhv 
ds ds ds 2 ds 2 ~ 


where $ denotes the angle of rotation of the cross section. 
Thus, the assumption of a linear variation of strain across the 
cross section yields 

dd> 

*ii + ( l6 ) 

Employing the principal of virtual work followed by integra- 
tion by parts and subsequent algebraic manipulation yields the 
equilibrium equations 

TV m -F^l + [a cos 0(0 + w(5) - w(/)]j (17a) 

M = F{a cos <j>{l ) .+ w(s) - vv(/)] ( 17b) 

where TV and A/, the force and moment resultants, respec- 
tively, are defined as 


N = 1 

| 

1 

- 

(18a) 

M = 1 

| VOudA 


(18b) 


A 


Based on an additive decomposition for the total strain. 

where e e and t c represent the elastic and creep 
strain components, respectively, yields, after substitution into 
Eqs. (17) and (18), " 


EAe 0 = -fjl +~\a cos <*>(/) + h(j) - »•(/)] j + N c (19a) 
d(t) 

El — = F [a co s0(/) + w(s ) - w(/)] + M c (19b) 
where the “pseudoresultants” TV C and M c are defined by 


N, = | Et c dA 

— (20a) 

A 


M c ~ | r]Et c dA 

(20b) 


A 
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Numerical solution lor Eqs. ( 19) arc computed using a 
modified Nevvton-tvpe method suggested by Thurston. To 
illustrate this method, consider a nonlinear differential term of 
the form dn m dic" where di< and die are differentials of the 
functions u and »• and m and n represent integer exponents 
Assutnine that close trial solutions if and 11 arc available, which 
differ from the exact solution by the small quantities Aw and A u 
so that u = ii 4- An and tf = w + Ate. then 


du'"d»'" = du'"dtf" + nidu'" 1 dit"d( An) 
+ iidii" l du ,n " 'd(Au) +/(ii,u)0[Aii.Au] 


( 21 ) 


where /( ) denotes a nonlinear function of u and tv and 0[ ] 
indicates terms of order An Air and higher. If the trial solution 
is indeed close to the true solution, then the corrections A u and 
An- will be small. Consequently, the quadratic and higher-order 
terms in the corrections will be negligible in comparison to the 
linear terms and therefore may be neglected. Thus, the left- 
hand side of Eq. (21) may be closely approximated by the 
linearized form consisting of just the first three terms on the 

right-hand side. . 

With this procedure, the original nonlinear differential equa- 
tion is approximated by a linearized form. Employing standard 
finite-difference formulas, the linearized form is then converted 
into a system of algebraic equations where the unknowns are 
the corrections to the trial solution at the nodes of the finite- 
difference mesh. These relations are solved for these correc- 
tions, the trial solution is adjusted, and the process repeated 
until convergence is obtained. 

Application of this procedure to the geometry of the sample 
problem, e.g., yields the finite-difference expressions 

EI 8 li A „ f _ , - 2E/ d -p A u, + A u l+ , 

ds - os os 

-Ei(a+ 2AS ^ A*v,_ , + (2 El - As 2 T) Aw, 
d 2 ii- 


- Ely 1 - 2As-^jJ Aw /+ ! + As 2 PAw n 
-f 2AsPa tan0(Aw n _ ! - w H + i)_ 


d 4>i 
ds 


= &s 2 M] +M C .-EI 


(22a) 


and 


M 


■ - 2Av£'.4 p.., - 2 AC ~ An, 


-«-(f 

-[2AsE A ^ + M‘(2As^ + l) 

+ 2A sP tan^jAw,^ j 4- As 2 P + Avv,- 


+ 2AsEA )Au / + 


dip, 

ds - 


As 2 PAw n 


+ 2A sP tan fa |Aw / + , 

-2 —2 AsPa tan<ji„ (A tv, _ , - A vv n + ,) 

3s 

= As 2 [ - P cos 4>i 4- N c . — £e 0 - ] 


where 


(22b) 

(23) 


In these equations, the subscript / is used to represent interior 
nodes, and the subscript n is employed to indicate the node at 
the loaded end of the beam-column. Values obtained from the 
trial solution have been denoted by the placement of a tilde 
over the applicable term. 

Numerical solution of Eqs. (22) requires evaluation of the 
“pseudoresultants” N c and M e at each interior point of the 
finite-difference mesh. This is accomplished by evaluating the 
accumulated creep strain at a select number of points across 
the cross section at each of the axial nodes. A three-point 
Newton-Cotes quadrature formula is then employed repeti- 
tively to approximate the area integrals. For the sample calcu- 
lations reported herein, evaluation of the accumulated creep 
strain at each of these points is accomplished through the use 
of a fourth-order Runge-Kutta integration routine to integrate 
the constitutive law. 


Example Problem 

The specific example considered is that of a 30.5 cm (12 in.) 
long beam-column. For simplicity, a square cross section ot 
dimension 12.7 mm (0.5 in.) has been assumed. It is also as- 
sumed that the beam-column is fabricated from a material that 
can be modeled as a three-parameter viscoelastic solid. The 
creep compliance for this model, illustrated in Fig. 2, is given 

by 

J(t)/J(0) = 1 +[£ 2 /£ 1 ]*~' /to ( 24 ) 

where t 0 = vJE { . For the sample computations, the numerical 
values for the parameters have been selected so that t 0 == 1. 
Thus, integer values for time t are equal to multiples of the time 
constant of the material. The elastic modulus of the material, 
E 2 , is assumed to be 196 GPa (28.5 x 10 6 psi), the room tem- 
perature modulus. of Hastelloy X. 

A five-point grid in the transverse direction is used in com- 
puting the “pseudoresultants” in the differential formulation. 
The points are equidistantly spaced, with the- first and last 
located at the extreme fibers and the central point positioned 
on the centroidal axis. 

Since the governing equations of both the differential and 
integral formulations are only satisfied at a discrete number of 
points over the length of the beam-column, the first question 
addressed is the sensitivity of the results to the number of 
points used in the approximation. Table 1, for example, com- 
pares initial elastic deflections determined using the differential 
formulation as the number of approximating points is doubled 
from 10 to 20 and then doubled again to 40. Table 2 provides 
a similar comparison for the integral formulation. 

It can be seen that there is very little change in the computed 
transverse deflection as the number of approximating points is 
increased. In both cases, the initial elastic solution for the 10- 
element model is within 1.0% of the 40-element model results. 
Additionally, the relative magnitude of the errors between the 
10- and 40- as well as the 20- and 40-eiement models of the 
differential formulation are very similar to those exhibited by 
the equivalent comparison of integral solution models. These 
specific results, of course, apply for an eccentricity ratio of 0.05 
and an applied load of P/P e ~ 0.75, where the Euler load P e is 
based on a perfect geometry and use of the instantaneous com- 
pliance of the material J( 0). However, they, like other results 
reported herein, illustrate the general trends observed at other 
load levels and load eccentricities. 

Thus, both formulations exhibit a similar low sensitivity to 
the number of elements used in the analysis. Concurrently, the 
results also indicate that a 10-element model can be used with 
either method without generating significant errors in the anal- 
ysis. It should be noted that all of the differential formulation 
results presented in the first table are based on the use of an 
exact expression for evaluating the angle of rotation. The influ- 
ence of employing an approximate formula lor calculating the 
angle of rotation is discussed in a later section. Also note that 


M* = P (a cos0„ 4. w, - vv„) 
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at litis load and eccentricity, the initial end rotation of the 
beam-column exceeds 17 deg. Of course, smaller rotations are 
exhibited at the lower loads and lower eccentricities. Higher 
loads and larger eccentricities, conversely, produce greater 
rotations. 

A direct comparison between the results generated by the 
two formulations is provided in Tables 3 and 4. Again, the 
comparison is based on the 0.05 eccentricity ratio, which pro- 
duces reasonably large angles of rotation. Angles of rotation of 
the cross section, determined from the integral technique, are 
included in the tabulated data. 

It should be noted that the differential solution results pre- 
sented in Table 3 are based on the use of an exact expression 
for evaluation of sin</>. Differential solution results obtained 
using an approximate expression for sin0 are provided in Table 
4. The common approximation <j> & sin" *( —dw/ds) is em- 
ployed to calculate the angle of rotation for this second set of 
results. Except for this particular difference, these two differen- 
tial formulations are otherwise completely identical. 

These results indicate that little difference exists between the 
initial deformation predicted using the integral formulation 
and that predicted by either differential solution. The differ- 
ences between the integral and differential methods are typi- 
cally an order of magnitude lower than the differences observed 
for either technique when the number of elements was quad- 
rupled. A potentially high-order effect may be indicated by the 
relative increase in differences at the highest loading examined. 
However, despite this increase, the magnitude of the differences 
is still so small as to be completely inconsequential with regard 
to engineering computations. 

The data also indicate that no significant differences in the 
differential formulation (predicted transverse deflections) oc- 
cur as a result of using an approximate formula for evaluation 
of sin<£. Even for an end rotation angle of 17 deg, the exact and 
approximate results differ only in the third or fourth decimal 
place. It should be noted that this high level of correlation 
continues to exist for even greater angles of rotation. 

The probable reason for this high correlation is that, under 
compressive loading, the derivative of the axial displacement is 
negative. With reference to Eq. (14), this implies that the cen- 
troidal axis strain is numerically equal to the difference be- 
tween the two components since the squared term (slope of the 
transverse deflection) is always positive. Thus, the magnitude 
of the centroidai axis strain must be less than the magnitude of 
either of its components. Because the difference between the 
ex act and approximate expressions for si ruf) is related to the 
y / 1 + 2e 0 in the denominator of the exact expression, reducing 
the magnitude of the centroidai strain must inherently improve 
the accuracy of an approximation where this term is neglected.- 

This is best illustrated by the data of Table 5. Here, the 
integral solution is compared to an approximate differential 
solution in which the effect of the centroidai axis strain terms 
are suppressed. This suppression is accomplished by eliminat- 
ing all the ( dw/ds ) 2 terms from the governing, equations. In 
addition, the EA modulus-area product is artificially increased 
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Fig. 2 Ideal three-element “limited” creep model. 


through multiplication by a factor of 1000. This second change 
reduces the magnitude of the axial deflections u by approxi- 
mately the same factor. The overall intent of this effort was to 
create a differential model that would simulate the axial “inex- 


Tablelf Differential elastic solution vs number of nodes for 
PjP t =* 0.75 and afl = 0.05 


Transverse deflection (cm) for various numbers of elements 


Number of elements 


s—s/t 

40 

20 

% diff 1 

0.0 

0.000000 

0.000000 

__ 

0.1 

0.062669 

0.062751 

0.13 

0.2 

0.249400 

0.249700 

0.12 

0.3 

0.556448 

0.557103 

0.12 

0.4 

0.977725 

0.978855 

0.12 

0.5 

1.504980 

1.506695 

0.11 

0.6 

2.128050 

2.130438 

0.11 

0.7 

2.835125 

2.838262 

0.11 

0.8 

3.613038 

3.616978 

0.11 

0.9 

4.447532 

4.452313 

0.11 

1.0 

5.323528 

5.329169 

0.11 

s — sfl 

40 

10 

% diff 

0.0 

0.000000 

0.000000 



0.1 

0.062669 

0.063015 

0.55 

0.2 

0.249400 

0.250858 

0.59 

0.3 

0.556448 

0.559534 

0.56 

0.4 

0.977725 

0.983259 

0.57 

0.5 

1.504980 

1.513218 

0.55 

0.6 

2.128050 

2.139785 

0.55 

0.7 

2.835125 

2.850330 

0.54 

0.8 

• 3.613038 

3.632429 

0.54 

* 0.9 

4.447532 

4.470819 

0.52 

1.0 . 

5.323528 

5.351315 

0.52 


•% differences are with respect to 40-element solution. 


Table 2 Integral elastic solution vs number of nodes for 
P/P t *0.75 and a// = 0.05 


Transverse deflection (cm) for various numbers of elements 


Number of elements 


s — sfl 

40 

20 

% diff 1 

0.0 

0.000000 

0.000000 

' 

0.1 

0.062525 

0.062616 

0.15 

0.2 

0.248836 

0.249238 

0.16 

0.3 

0.555208 

0.555981 

0.14 

0.4 

0.975561 

0.977115 

0.16 

0.5 

1.501648 

■ 1.503906 

01.5 

0.6 

2.123316 

2.126658 

" 0.16 

0.7 

2.828762 

2.833050 

0.15 

0.8 

3.604837 

3.610412' 

0.16 

0.9 

4.437309 

4.443969 

0.15 

1.0 

5.311148 

5.319194 

0.15 

J - s/I 

40 

10 

% diff 

0.0 

0.000000 

0.000000 



0.1 

0.062525 

0.062944 

0.67 

0.2 

0.248836 

0.250833 

0.80 

“ 0.3 

0.555208 

0.5591 12 

0.70 

0.4 

0.975561 

0.981829 

0.64 

0.5 

1.501648 

1. 5 1 2425 

0.72 

0.6 

2.123316 

2.137987 

_0.69 

0.7 

2.828762 

2.847177 

0.65 

_ 0.8 

3.604837 

3.629746 

0.69 

0.9 

4.437309 

4.466392 

0.66 

1.0 

5.311148 

5.348440 

0.70 


*% differences are with respect lo 40-e!emenl solution. 


38 



Tabic 3 Comparison* of integral and dilTcrential b elastic solotions for 
various loads with ajl = 0.05 


Transverse deflection (cm) from various solutions 


=* Sil 

Integral 

Differential 

% difT 

Angle, 

cleg 

0.0 

0.000000 

PIP, = 0.25 
0.000000 

_ 

0.00 

0.1 

0.006642 

0.006645 

0.04 

0.25 

0.2 

0.026543 

0.026538 

-0.02 

0.50 

0.3 

0.059548 

0.059550 

0.00 

0.74 

0.4 

0.105461 

0.105489 

0.03 

0.98 

• 0.5 

0.164059 

0.164056 

0.00 

1.22 

0.6 

0.234887 

0.234907 

0.0 1 

1.44 

0.7 

0.317525 

0.317579 

0.02 

1.66 

0.8 

0.411579 ■ 

0.411592 

0.00 

1.87 

0.9 

0.5 16270 

0.516329 

0.0 1 

2.07 

1.0 

0.631223 

0.631182 

-0.01 

2.25 

0.0 

0.000000 

PIP , = 0.50 
0.000000 


0.00 

0.1 

0.021039 

0.021044 

0.02 

0.79 

0.2 

0.083962 

* 0.083932 

-0.04 

1.57 

0.3 

0.187808 

0.187828 

0.01 

2.33 

0.4 

0.331320 

0.331511 

0.06 

3.07 

0.5 

0.513103 

0.513077 

0.00 

3.76 

0.6 

0.730283 

0.730400 

0.02 

4.40 

0.7 

0.980300 

0.980618 

0.03 

5.01 

0.8 ‘ 

1.260747 

1.260810 

0.01 

5.52 

0.9 

1.566974 

1.567304 

0.02 

6.01 

1.0 

1.896801 

t. 896542 

-0.01 

6.38 

0.0 

0.000000 

P/P e = 0.75 
0.000000 


0.00 

0.1 

0.062944 

0.063015 

0.11 

2.37 

0.2 

0.250833 

0.250858 

' 0.01 

4.68 

0.3 

0.5591 12 

0.559534 

0.08 

6.91 

0.4 

0.981829 

0.983259 

0.15 

9.03 

0.5 

1.512425 

1.513218 

0.05 

10.97 

0.6 

2.137987 

2.139785 

0.08 

12.68 

0.7 

2.847177 

2.850330 

0.1 1 

14.23 

0.8 

3.629746 

3.632429 

0*07 

15.41 

0.9 

4.466392 

4.470819 

0.10 

16.44 

1.0 

5.348440 

5J51315 

0.05 

17.04 


•Comparisons based on results from 10-element models, differential solution 
employing an exact sin<^ formula. c % differences are with respect to integral 


solution. 

tensionality” of the integral model. These changes did produce 
a differential model with an effectively inextensional centroidal 
axis. It was anticipated that this would further improve the 
correlation between the differential and integral results. Unfor- 
tunately-, such was not the case. 

When the angle of rotation is very small, such as one that 
results from a low level of loading and minimal eccentricity, all 
formulations provide virtually identical predictions. Increases 
in the angle of rotation, however, due to increases in loading or 
eccentricity or both, cause the modified differential predictions 
to diverge from those of the others. This divergence between 
results increased with both load magnitude and eccentricity. 

This behavior is attributed to the manner in which the 
modified numerical model handles the end deflection of the 
beam-column. In the modified model, the end of the beam- 
column effectively moves only in the vertical direction (see Fig. 
1). The standard differential model as well as the integral 
model, however, include influences generated when the end can 
move both vertically and horizontally. Thus, for any given 
vertical deflection, the horizontal movement that occurs in the 
integral and unmodified differential models acts to increase the 
angle of rotation. This, in turn, reduces the magnitude of the 
applied moment [see Eq' (23)]. Therefore, at any particular 
given vertical deflection, the moment loading in the standard 
formulation model is lower than that in the modified version. 
Effectively, the moment decreases a greater amount in the un- 


Tabie 4 Comparison- of integral and ainerennar elastic solutions lor 
various loads with aji ~ 0.5 

Transverse deflection (cm) from various solutions 


Angle, 


s — sjl 

Integral 

Differential 

% difF 

deg 

0.0 

0.000000 

PIP, = 0.25 
0.000000 

_ 

0.00 

0.1 

0.006642 

0.006645 

0.00 

0.25 

0.2 

0.026543 

0.026538 

-0.02 

0.50 

0.3 

0.059548 

0.059550 

0.00 

0.74 

0.4 

0.105461 

0.105489 

0.03 

0.98 

0.5 

0.164059 

0.164056 

0.00 

1.22 

0.6 

0.234887 

0.234907 

0.01 

1.44 

0.7 

0.317525 

0.317579 

0.02 

1.66 

0.8 

0.411579 

0.411592 

0.00 

1.87 

0.9 

0.516270 

0.516329 

0.01 

2.07 

1.0 

0.631223 

0.631182 

-0.01 

2.25 

0.0 

0.000000 

PIP, = 0.50 
0.000000 


0.00 

0.1 

0.021039 

0.021046 

0.04 • 

0.79 

0.2 > 

0.083962 

0.083932 

-0.04 

1.57 

0.3 

0.187808 

0.187830 

0.01 

2.33 

0.4 

0.331320 

0.331511 

0.06 

3.07 

0.5 

0.513103 

0.513080 

0.00 

3.76 

0.6 

0.730283 

0.730402 

0.02 

4.40 

0.7 

0.980300 

0.098620 

0.03 

5.01 

0.8 

1.260747 

1.260813 

0.01 

5.52 

0.9 

1.566974 

1.567307 

0.02 

6.01 

1.0 

1.896801 

1.896547 

-0.01 

6.38 

0.0 

0.000000 

P/P € = 0.75 
0.000000 


0.00 

0.1 

0.062944 

0.063017 

0.12 

2.37 

0.2 

0.250833 

0.250868 

0.01 

4.68 

0.3 

0.559112 

0.559559 

0.08 

6.91 

0.4 

0,981829 

0.983305 

0.15 

9.03 

0.5 

1.512425 

1.513286 

0.06 

- 10.97 

0.6 

2.137987 

2.139881 

0.09 

12.68 . 

0,7 

2.847177 

2.850459 

0.12 

14.23 

0.8 

3.629746 

3.632594 

0.08 

15.41 

0.9 

4.466392 

4.471020 

0.10 

16.44 

1.0 

5.348440 

5.351554 

0.06 

17.04 


•Comparisons based on results from 10-eicment models, differential solution 
employing an approximate sin <f> formula. c % differences are with respect to 


integral solution. 


modified model than it does in the modified model at equiva- 
lent amounts of transverse deflection. This, in turn, implies that 
the beam-column of the unmodified model would not deflect as 
much as the one of the modified model would. 

The implication is that the numerical modeling of the influ- 
ence of deformation on loading is an important factor. This 
conclusion is consistent with the observations made concerning 
how the structure of the integral equation, Eq. (5), changes as 
a result of the interaction between loading and deflection. Ad- 
ditionally, it should be noted that neglecting centroidal axis 
strain in'the differential technique does not necessarily provide 
an effect equivalent to the assumption of inextensionality in the 
integral techniques. This is attributed to the fact that the elim- 
ination of the centroidal axis strain in the differential formula- 
tion can be accomplished only at the expense of reducing the 
actual coupling between deformation and loading. 

For both methods, the initial elastic deffection of the beam- 
column provides the initial condition for the viscoelastic defor- 
mation. Consequently, any differences in the initial elastic 
responses predicted by the two methods will only be accentu- 
ated during the subsequent period of time-dependent behavior. 
The comparisons discussed demonstrate that the two methods 
provide virtually identical predictions of initial clastic deforma- 
tion. Table 6 provides a typical comparison for the integral and 
the differential (exact sin <f>) method predictions for the vis- 
coelastic deflection of the beam-column over a period of two 
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Table 5 Comparison* of integral and modified differential 6 elastic 
solutions for various loads with a/l = 0.05 


Transverse deflection (cm) from various solution 


Angle, 


s = si! 

Integral 

Differential 

% difF 

deg 

0.0 

0.000000 

P/P, = 0.25 

o.oooooo 


0.00 

0.1 

0.006642 

0.006645 

0.00 

0.25 

0.2 

0.026543 

0.026533 

-0.04 

0.50 

0.3 

0.059548 

0.05954S 

0.00 

0.74 

0.4 

0.105461 

0.105479 

0.02 

. 0.98 

0.5 

0.164059 

0.164048 

-0.01 

1.22 

0.6 

0.234887 

0.234894 

0.00 

1.44 

0.7 

0.317525 

0.317576 

0.02 . 

1.66 

0.8 

0.411579 

0.411 589 

0.00 

1.87 

0.9 

0.516270 

0.516349 

0.02 

2.07 

1.0 

0.631223 

0.631210 

0.00 

2.25 

0.0 

0.000000 

P/P, = 0.50 
0.000000 


0.00 

0.1 

0.021039 

0.021067 

0.13 

0.79 

0.2 

0.083962 

0.084005 

0.05 

1.57 

0.3 

0.187808 

0.188041 

0.12 

2.33 

0.4 

0.331320 

0.331889 

0.17 

3.07 

0.5 

0.513103 

0.513776 

0.13 

3.76 

0.6 

0.730283 

0.731457 

0.16 

4.40 

0.7 

0.980300 

0.982246 

0.20 

5.01 

0.8 

1.260747 

1.263051 

0.18 

5.52 

0.9 

1.566974 

1.570406 

0.22 

6.01 

1.0 

1.896801 

1.900517 

0.20 

6.38 

0.0 

0.000000 

PIP, = 0.75 
0.000000 


0.00 

0.1 

0.062944 

0.064879 

3.07 

2.37 

0.2 

0.250833 

0.258313 

2.98 

4.68 

0.3 

0.559112 

0.576725 

3.15 

6.91 

0.4 

0.981829 

0.014219 

3.30 

9.03 

0.5 

1.512425 

1.562702 

3.32 

10.97 

0.6 

2.137987 

2.212025 

3.46 

12.68 

0.7 

2.847177 

2.950167 

3.62 

14.23 

0.8 

3.629746 

3.763475 

3.68 

15.41 

0.9 

4.466392 

4.636892 

3.82 

16.44 

TO 

5.348440 

5.554259 

3.85 

17.04 


•Comparisons based on results from 10-element models, influences of centroidai 
axis strain terms suppressed. c % differences are with respect to integral solution. 


material time constants. The viscoelastic model employed for 
these computations is the three-parameter “limited” creep ma- 
terial illustrated in Fig. 2. As noted previously , the material 
parameters were selected so that the material time constant t 0 
equals unity. The load and eccentricity ratios for this particular 
set of results were P/P e = 0.50 and a/l = 0.05, respectively. 

As demonstrated by this data, the high correlation between 
the integral and differential method predictions for the initial 
elastic deflections carries over directly to the viscoelastic analy- 
sis. The time-dependent deflection predicted by one technique 
is virtually indistinguishable from that predicted by the other. 
This indicates that the differential formulation methodology 
employed to account for the influence of viscoelastic strain 
provides the equivalent effect as the hereditary integral compo- 
nent of the integral formulation. As such, this lends high confi- 
dence to the differential solution methodology. 

It should be noted that these particular numerical results are 
typical of other results obtained for higher, as well as lower, 
loads and eccentricities. Generally, the correlation between the 
solutions was not influenced by the magnitude of the loading or 
the amount of load eccentricity. 

A minor, high-order-type influence was, however, noted. As 
the angle of rotation became very large, on the order of 45 deg, 
a small but distinct divergence in predicted deflections was 
observed. Typically, the rate of increase in deflection predicted 
by the differential formulation would begin to slightly exceed 


Table 6 Comparison* of integral and differential 6 viscoelastic solutions 
to two time constants; P/P, =0.50 and ajl =0.05 

Transverse deflection (cm) for various solutions 


Angle, 


s * sfl 

Integral 

Differential 

% difF 

deg 

0.0 

0.000000 

at time = 0 
0.000000 


0.00 

0.1 

0.021039 

0.021044 

0.02 

0.79 

0.2 

0.083962 

0.083932 

-0.04 

1.57 

0.3 

0.187808 

0.187828 

0.01 

2.33 

0.4 

0.331320 

0.331511 

0.06 

3.07 

0.5 

0. 513103 

0.513077 

0.02 

3.76 

0.6 

0.730283 

0.730400 

0.02 

4.40 

0.7 

0.980300 

0.980618 

0.03 

5.01 

0.8 

1.260747 

1.260810 

0.0 1 

5.52 

0.9 

1.566974 

1.567304 

0.02 

6.01 

1.0 

1.896801 

1.896542 

-0.01 

6.38 

0.0 

0.000000 

at time = r 0 
0.000000 



0.00 

0.1 

0.038311 

0.038326 

0.04 

1.44 

0.2 

0.152784 

0.152725 

-0.04 

2.86 

0.3 

0.341170 

0.341234 

0.02 

4.22 

0.4 

0.600532 

0.600994 

0.08 

5.54 

0.5 

0.927608 

0.927590 

0.00 

6.76 

0.6 

1.315809 

1.316129 

0.02 

7.86 

0.7 - 

1.759290 

1.760083 “ 

0.05 

8.87 

0.8 

2.252647 

2.252896 

0.0 1 

9.70 

0.9 

2.785532 

2.786395 

0.03 

10.44 

1.0 

3.353191 

3.352810 

- 0.01 

10.95 

0.0 

0.000000 

at time = 2t 0 
0.000000 


0.00 

0.1 

0.048583 

0.048611 

0.06 

1.83 

0.2 

0.193685 

0.193629 

—0.03 

3.62 

0.3 

0.432178 

0.432313 

0.03 

5.34 

0.4 

0.759968 

0.760689 

0.09 

7.00 

0.5 

1.172517 

1.172627 

’ 0.01 

8.52 

0.6 

• 1.660764 

1.661396 

0.04 

9.88 

0.7 

2.216699 

2.218045 

0.06 

11.13 

0.8 

2.832971 

2.833680 

0.03 

12.12 

0.9 

3.495617 

3.497252 

0.05 

12.99 

1.0 

4.198267 

4.198371 

0.00 

13.55 


“Comparisons based on results from 10-element models, differential solution 
employing an exact sin<£ formula. c % differences are with respect to integral 
solution. 

that predicted by the integral method. Normally, this could be 
observed as time approached two material time constant for 
the highest loads (P/P e approaching unity) and with extremely 
large eccentricities. Under some conditions, it also could be 
observed as time exceeded four to five material time constants. 

These observed differences between the two sets of predic- 
tions were still very small. Generally, they were on the order of 
0.10%. Thus, from a practical viewpoint, they are totally negli- 
gible with respect to normal accuracy requirements for engi- 
neering computation. It is mentioned here only to indicate that, 
under conditions such as these, the accumulation of numerical 
errors may begin to influence the results. A comparison be- 
tween the integral and the approximate s\ruj> differential solu- 
tion results was not included because the differences between 
the exact and approximate differential solution results are 
again so small as to be negligible. 

The final item meriting discussion is the length of the time 
increment used in each of the formulations. Unlike the prior 
results, some differences do exist between the maximum allow- 
able time-step increments for the integral and differential for- 
mulations. Additionally, the allowable time-step increment for 
the integral formulation exhibits a higher dependence on the 
actual angle of rotation than does the differential formulation. 

In general, a relatively small time step increment must be 
used with the integral solution methodology. For example, the 
results previously presented typically employed a 0.01 -time- 
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step increment. As the length of this time step is increased, the 
accuracy of the solution decreases and tends to underprcdict 
the deflection. This convergence "from below" is not surprising 
since ihc convolution integral is approximated as the sum of a 
finite number of terms. 

In contrast, much larger time-step increments were used with 
the differential formulation. This is principally attributed to the 
hich accuracy provided by the Runge-Kutta integration rou- 
tine. Most of the results provided were developed using a 0.10 
time step. The use of even larger time steps was also examined. 

It was found that time increments four to five times greater 
than 0. 10 could be employed without significant changes in the 
calculated results. Additionally, the allowable length of this 
time-step increment tended to be rather insensitive to the angle 
of rotation. The allowable time step for the integral formula- 
tion. on the other hand, exhibited a high level of sensitivity to 
the angle of rotation. Larger angles of rotation required signifi- 
cantly shorter time steps for accurate results to be obtained. 

These factors combine in a rather interesting manner with 
regard to which method of analysis is computationally more 
efficient. Typically, for the analysis of short periods of vis- 
coelastic deformation, the integral solution method was two to 
three times faster than the differential method. This is at- 
tributed to two factors. The first is the comparitively slow 
Runge-Kutta integration procedure used in the differential for- 
mulation. For a short period of viscoelastic deformation, the 
calculation of the convolution integral of the integral formula- 
tion, requiring simple summation of a limited number of terms, 
can be performed much more rapidly. 

The second factor is that the fixed-point iteration scheme of 
the integral formulation, although requiring more iterations 
than the Newton method, is also performed more rapidly since 
it is simply an algebraic operation. The Newton method, in 
contrast, requires inversion of the matrix, premuitiplying the 
vector of trial function corrections, and then numerical evalua- 
tion through solution of the system of equations. Even for just 
a 10-element beam, this process is slow in comparison to the 
fixed-point iteration. . 

However, as the length of the period of viscoelastic deforma- 
tion increases, this relative speed relationship reverses. Eventu- 
ally, the differential formulation begins to generate solutions 
more rapidly than the integral method. In the example problem 
previously described, this generally occurred approximately be- 
tween the second and third time constants. The reason for this 
change is directly related to the computation of the convolu- 
tion integral of the integral technique. As time increases, the 
number of terms in the.summation increases linearly. This, in 
turn, increases the number of algebraic operations that must be 
performed and therefore linearly increases the time need for 
each complete computation.. In contrast, the speed of the 
Runge-Kutta integration routine is virtually independent of 
time. Thus, the continually increasing computational effort re- 
quired in the integral technique eventually exceeds that need 
for the differential technique. This reverses the relative speed 
relationship. 

Conclusions 

Based on the results reported herein and elsewhere, 12 it is 
concluded that the differential formulation procedure pre- 


sented can be employed for the analysis of quasistatic non- 
linear one-dimensional viscoelastic problems. This conclusion 
is based directly on the high level of correlation between results 
developed using this formulation technique to those obtained 
with the previously published integral method for solution o 
such problems. Additionally, it is observed that both of these 
methods exhibit exceptionally similar accuracy characteristics 
with resard to the number of elements employed in the approx- 
imation. For both, a relatively low number of elements can be 
used without engendering any significant errors. 
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HEREDITARY V1SCO- ELASTIC- PLASTIC CONSTITUTIVE LAW 
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ABSTRACT 

An analytic study of planar beams and arches subjected to signifi- 
cant thermal cycling from ambient temperatures up to 800 "C is presented. 
In the study, a recently developed unified nonlinear hereditary type of vis- 
eoel&stoplastic constitutive law is employed to characterise the time- and 
temperature-dependent properties of a typical aerospace alloy, Ha&telloy X. 

The results from this work demonstrate that a strong interaction 
exists between the backstress variable of this particular constitutive law 
.and the time-dependent stress distribution produced by the geometry of 
the deformation. Effectively, this interaction tends to control, in, a highly 
nonlinear manner, the creep-ratchetting response of the beam and the arch. 
An unexpected consequence of this is that temperature gradients in the 
thickness direction, a factor normally neglected in most studies, tends to 
exert an important influence on the response during thermal cycling. 


NOMENCLATURE 


A 

b 

d 

E 

E 0t Ei , E j 


g«» G« 
9*0* G*0 


h 

k 

K 


K U K* 

m thru nr 


m, n 


load eccentricity 

cross sectional area of beam or arch 
width of beam or arch 
inelastic strain tensor 

deformation vector for points on the: centroidal axis 
Young’s Modulus 

sero, first and, second moments.. of the elastic modulus 
across a cross section, respectively 

base vectors for undeformed and deformed configurations 
metric components of the undefornied and oeformed 
configurations 
depth of beam or arch 
*+ 

drag stress 

constitutive law constants 
constitutive law constants 
constitutive law exponents 


M, N 
M e , N e . 
M#, N# 
p 

P, V 
Q(t)> T(t) 
r, R 
** V 

t, n, k 
T, N, K 
i 

u, u> 

Of 

e # 

Uj 

© 

K 

A, m 

<Tij 


moment and force resultants 

creep strain moment and force pseudo-resultants 
thermal strain moment and force pseudo-resultants 
pressure load 

axial and transverse force resultants 
constitutive law functions 

position vectors in undefornied and deformed configurations 

coordinates along length and depth directions 

components of the deviator stress tensor 

triad of unit vectors for the undeformed configuration 

triad of unit vectors for the deformed configuration 

time 

axial and transverse displacement of the centroidal axis 

coefficient of thermal expansion 

centroidal axis strain 

strain tensor component 

change in temperature 

initial curvature of the arch 

Lame constants 

stress tensor component 

angle of rotation of crosa section 

backstress 


INTRODUCTION 

It is well known that metal alloys can undergo transitions in behavior 
as temperature increases. Commonly, for loading substantially below yield, 
the elastic response observed at room temperature generally gives way to a 
time-dependent viscoelastic response at somewhat elevated temperatures. 
Further increases in temperature, however, introduces the potential for sud- 
den “rapid” or plastic type deformation. Such transitions can significantly 
shorten the useful life of the structural, element and generate the possibil- 
ity of a sudden unanticipated failure. Consequently, for many years the 
aerospace and nuclear power industries, where elevated temperature oper- 
ating, environments abound, have had a continuing interest in predicting 
the behavior of metallic structural elements subjected to such conditions. 

Early investigators [1-3] generally focused their attention on the be- 
havior of structural components subjected to conditions of constant load at 


42 



constant uniform elevated temperatures. Many of these studies employed 
simplified analyses to improve mathematical tractability. Additionally, ’‘ex- 
perimentally based equation of state” type constitutive laws were often used 
to express the nonlinear elevated temperature time-dependent behavior of 
the material. A summary of many of methods developed and key findings 
obtained is provided by Hoff [4]. 

These efforts answered many of the questions regarding elevated tem- 
perature creep buckling. However, they were not able to satisfactorily 
describe the creep ratchetting behavior resulting from thermal cycling at 
elevated temperatures. Consequently, it was not until Miller [5] and Ed- 
munds and Beer [6], both of whom considered nuclear pressure vessels, that 
this particular form of behavior was specifically addressed. Subsequently, 
Bree [7,8] investigated basic factors which determine when this type of re- 
sponse could occur. These studies led to experimental investigations and 
other analytic studies to explore various aspects of the phenomena. The 
works of Conway et. al. [9], Corum [10,11] and Mukhejeree, Kumar and 
Chang [12] provide a representative sampling of these efforts. 

However, the greatest concentration of effort has been directed to- 
ward improving the capability to predict the elevated temperature behav- 
ior of metals. The contributions of Hart [13,14], Pointer and Leckie [15], 
Pugh [16-18], Krempl [19], and Walker and Krempl [20], to name of few, . 
provide a dramatic illustration of the intensity of these efforts to develop ad- 
vanced constitutive models. Yet, as pointed out by Corum and Sartory [21], 
an equation of state approach to constitutive modeling is still generally used 
in design situations. However, as Pugh [18] has noted, the newer types of 
unified constitutive laws, where inelastic strain is not divided into distinct 
creep and plasticity components, can provide an alternate approach. 

Consequently, one of the aims of this study is to examine the use of a 
typical unified constitutive model in an analysis of the behavior of structural 
elements subject to thermal cycling from an elevated ambient temperature. 
The specific law which was selected is one developed by Walker [22] to 
model the time- and temperature-dependent behavior of Hastelloy X, an 
alloy routinely used in the aeropsace industry. 

The study results indicate that, with this particular constitutive law 
and material, an implicit interaction exists between the stress in the mem- 
ber and the backstress of the constitutive law. This interaction strongly 
influences the ultimate response. This result is significant for two reasons. 
First, because saturation of the backstress of this constitutive law can lead 
to plastic response, it indicates that the ultimate reliability of any predicted 
results rests strongly upon. the accuracy with which the backstress growth 
law parameters have been determined. 

Additionally, this aspect produces the rather interesting result that, 
due to the strong temperature dependence of the material constants of 
the law, temperature variations in the thickness direction greatly influence 
predicted response. This is significant for any analysis since the influences 
of such temperature variations are generally neglected. 

MATHEMATICAL FORMULATION 

To focus principally upon the interaction between the response of 
the structural element and the prediction of thermal dependence of theana- 
terial, the problem is formulated within the context of a simplified beam 
theory. Consequently, it is assumed that the beam or arch deforms in accor- 
dance with the Euler- Bernoulli hypotheses. As such, cross sectional planes 
normal to the centroidal axis in the undeformed geometry are assumed 
to remain plane and normal in the deformed state. Similarly, extensional 
straining in the thickness and depth directions are neglected. Thus, based 
on the geometry illustrated in Fig. 1, the position vectors r and R, where 

r = r # + rjn and R = r, 4* d + 0 ) 

are employed to locate a typical point on an arbitrary cross section in the 
undeformed and deformed configurations, respectively. Note that r) repre- 
sents the coordinate in the normal direction. Also, lower case and upper 


case symbols are employed to denote quantities referred to the undeformed 
and deformed configurations, respectively. 

Base vectors for the reference and current state, g a and G„> respec- 
tively, are defined by 


Consequently, the deformation vector which translates a point from the 
undeformed to deformed configurations, denoted as d, can be expressed as 


d + >?N = (u + T?sin(£)t -f (u> + r?cos<p)n ( 3 ) 

where u and w represent axial and transverse displacement functions for 
points on the centroidal axis, respectively. Substitution of Eqns. (l) and 
(3) into (2), followed by differentiation and subsequent employment of the 
Fernet-Serret formulae and the strain definition, 


9ae) i ff./> = 6«*S0. Go/) - G “ • G 0 (4) 

yields the strain expressions 

1 .. = i{(l + ^ + IC1U + n*COS^) , + (^-Kt‘-’?* Sin ^) , ~( 1 + VC ) 3 } ^ 


= i {(1 + + Kw + cos 4>) sin <f> + — KU — T)k sin <£) cos^i} (5fc) 

and 

7uu = i{sin J ^ + cos J <£- l} = 0 ( Sc ) 

In these, k - * + 6<f>/da , where k denotes the initial constant curvature of 
the arch and <f> represents the.angle of rotation of the cross section. 

From the Euler-Bernoulli hypotheses, the shear strains must vanish. 
Consequently, from this requirement Eqn. (5b) yields 


f dw v 

— du : 
(!+ +**,) 


Therefore, employing the definition 


, du /OUf 

1 + 2e* = (l + + kw ) + (jj - **) 


(6) 


(?) 



Figure 1. Geometry of Deformation. 
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and Eqn. (6) it is relatively easy to show that 


sin 4> = 


, dw 

'77 


») 


s/\ + 2«. 


and cos^ = 


9u 

1 + - r - + kw 

os 

y/i + 2<» 


( 8 ) 


Thus, expanding Eqn. (5a) and employing Eqns. (7) and (8) yields 


_ , 9<f> > j, /c\ 

7,« = + it* -K+Vl+2£»^}+7) + 2 Os* ds 

From the above it is clear that < a represents the strain induced alon^ the 
centroidal axis. Since, for a thin arch or beam, the last term of Eqn. (9) 
should be small in comparison to the others, it may be neglected. Similarly, 
additional simplifications may be obtained for the case where the centroidal 
axis strain is sufficiently small so that it may be neglected in comparison to 
one. Based on these assumptions, Eqn. (9) simplifies to the standard form 


7,. ««• + ») 


ds 


( 10 ) 


Equilibrium equations are obtained through application of the Prin- 
ciple of Virtual Work. Stress resultants, N and M, are defined such that 


N= // 


Note that 


and M = / <rr)dA 

(u) 

Ja 


d6v\ , e 
— - + kSv 2 

(12o) 

9s 


96v 2 . f 

(126) 


where and 6t denote the incremental changes in rotation and centroidal 
axis strain resulting from the deformed configuration axial and transverse 
displacement changes, Sv\ and 6v 2t respectively. 

Consequently, from the Principle of Virtual Work, 


/ ! 
-p'Sv 7 d$ (13) 

Expressing the first term on the right hand side in terms of an integral over 
the length and then combining that result with the work term yields 
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Therefore, employing Eqns. (12) and noting that terms mult.plymg the 
virtual displacements f» t , Sv 2 and 6<t> must vanish identically yields, after 
eliminating the shear resultant, the equilibrium equations 
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with the associated boundary conditions at s - 0, / 



or tfvi — 0 

or Sv 2 = 0 (1®) 

or &<f> = 0 

where N\ Q\ and M* denote the axial, shear and moment resultants ap- 
plied at the ends of the arch, respectively. i 

Expressions for the force and moment resultants, N and M, respec- 
tively, are obtained from the constitutive law. A unified hereditary visco- 
elastic-plastic law developed by Walker [22] to characterise the time and 
temperature dependence of Hastelloy X is employed. The selection of 
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Figure 2 External Forces and Moments Acting on the Arch. 

Walker’s functional theory, a highly generalized representation for a three 
parameter viscoelastic solid, was based on three considerations. F.rst 
from a strictly practical point of view, a substantial body of experimental 
work [22-24] had been performed to establish the temperature dependence 
of the constitutive law parameters for a wide range of temperatures. These 
efforts included a validation which examined the predictive capab.hty of the 
law through a comparison of analytical and experimental results for time 
variable thermo-mechanical load cycling of uniaxial specimens. 

A second factor which favored selection of the Walker law was that it 
is able to reproduce forms of classical behavior as limiting cases. For exam- 
ple, saturation of the drag stress produces an effect equivalent to isotropic 
hardening of a material. Similarly, saturation of the backstress produces an 
effect equivalent to kinematic hardening. Finally, the associated laws which 
govern the evolution of the state variables provide for the opportunity to 
include effects related to both dynamic and static thermal recovery. 

The final reason for selection of the Walker law is that it can be ex- 
pressed in both differential and integral formats. For this particular study, 
the differential format of the law was found to be the most convenient. 
However, it was hoped that the availability of an integral format would 
provide, at a future date, the opportunity to extend some prior work in- 
volving kinematic bounding of nonlinear integral formulations [25] to also 
include some form of constitutive bounding. 

The general integral form of Walkeris functional theory has been 
provided in Appendix A. That appendix also contains a derivation of the 
differential form from the integral format. It should be noted that, in the 
modeling of the elevated temperature behavior of Hastelloy X, additional 
simplifications in the form of the law were possible. These simplifications 
resulted from the fact that, for Hastelloy X, a number of material constants 
are xero over the entire temperature range. 

Due to these simplifications, the differential format of Walker’s func- 
tional theory, for one-dimensional loading of Hastelloy X, has the form 

it = sign(<r - flu) ( 17a ) 

nu = oi l + n J <.-(n„-n; l )(6-i^©) (n&) 

n t , = ^e (nc) 


K = Ki (lTd) 

and 

G = n s abs(c* c ) + n 6 (abs{nu}) m 1 (17e) 

Finally, since the time rate of change of temperature is relatively 
low in the sample problems, all terms where © appeared were assumed 
to be insignificant and therefore neglected. Note that this provides some 
minor simplifications to Eqn. (17b) and makes the reference backstress, 
flJ lP independent of the time rate of change in temperature. 

The law is based on an additive decompostion of the strain into 
elastic, thermal and inelastic components, and <«, respectively. Thus, 
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For the state 


of one-dimensional loading considered, this yie'd* 
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( 18 ) 


(19) 


Therefore, integrating this expression over the cross section and employing 
the force resultant definition, Eqn. (11), yields 

( 20 ) 

where B 0 and E, denote the zero and first integrals of the elastic modulus 
over the cross section and N, and N« represent “pseudo-resultant” type 
quantities defined as 


n jl 
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( 21 ) 


Note that the first integral of the elastic modulus over the cross section, Bi , 
does not necessarily vanish since the elastic modulus is a function of the 
temperature, which is not necessarily constant across a cross section An 
interseting aspect of this is that it induces a form of bending-stretching 
coupling {see also Eqn. (22)} similar to that of a laminated compos, te 
material. Finally, it should be noted that, although the quantities defined 
by Eqn. (21) have the units of a force resultant, these definitions are merely 
employed to simplify subsequent expressions. 

Multiplying Eqn. (19) by the coordinate, tj, integrating over the cross 
section and employing the moment resultant definition, Eqn. (12), yields 

M = + ( 22 ) 

' where E 2 denotes the second integral of the elastic modulus over the cross 
section and the “pseudo-moment” resultants are defined by 

dA and M c = f r)Ec e dA (23) 
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Consequently, substituting Eqns. (20) and (22) into Eqns. (15) yields the 
general governing relations 


and 


The basic methodology is to expand a nonlinear differential term 
such as dX m dY", where dX and dY represent differentials of the functions 
X and Y and m and n represent integer exponents, into the sum and prod- 
ucts of trial solutions X and Y and corrections of the form AX and AY. 
Consequently, substituting X = X -f AX and Y = Y + AY yields, 


dX m dY" = dX m dY n + max”*" 1 dY n d(AX) 

+ ndX m dY n ” l d(AY) + /(X, Y)0[AX, AY] 


( 26 ) 


(& + wife . (* + *&>£ = £<»•+ N '» + * (M - + “*> (!S) 

In the development of these relationships, it is assumed that the moments of 
the elastic modulus, Bo, El and E 3 , are constant with respect to the axial 
coordinate. This implies that these quantities are independent of stress or 
strain and that the temperature is constant along the length direction. 

Equations (24) and (25) represent a pair of interrelated spatially de- 
pendent nonlinear differential relationships which describe the deformation 
of the beam or arch. In their present form, they are stated in terms of the 
strain along the centroidal axis and the rotation of the cross section. These 
two quantities, in turn, are interrelated through the axial and transverse 
displacements of the centroidal axis. Substitution of the appropriate ex- 
pressions for the centroidal axis strain and cross section rotation ultimately 
results in a set of equations, one of third order and the other of fourth 
order, in terms of these displacement functions. 

Due to the significant nonlinearity of these equations, a numerical 
method of solution was selected. The particular method employed is an 
adaptation of Newton’s method for the solution of nonlinear algebraic equa- 
tion. [26]. The basic approach is an iterative procedure where a “close trial 
solution is directed towards the actual solution. Note that this method nei- 
ther guarantees convergence to « solution nor that a solution is unique. 


where /(X, Y)0(AX, AY) represents a nonlinear function of X and Y of 
second and higher order terms in AX and AY. Provided the trial solution 
is close to the true solution, these higher order terms should be small in 
comparison to the linear terms and thus may be neglected. Consequently, 
the nonlinear term may be closely approximated by the linear form 

dX m dY“ « dX m dY n + rndX m -* dY"d(AX) + ndX m dY"' 1 d( AY) (27) 

With this technique, the nonlinear differential equations are approxi- 
mated in terms of linear differential equations for the corrections to assumed 
trial functions. These coupled differential equations arc then converted to 
a set of coupled algebraic relations through the use of central difference 
formulae to approximate the derivative terms for the corrections to the as- 
sumed deflections. Appendix B provides, for example, the general finite 
difference expressions developed for the initially circular arch. 

A matrix iteration procedure is employed to refine an assumed trial 
solution. Each sucessive set of corrections is used to update the trial solu- 
tion until convergence is obtained. Tests for such convergence included the 
consideration of the magnitude of each set of corrections as well as the over- 
all accuracy for which each of the individual nodal equations was solved. 

In this regard it is noted that an equivalent degree of coupling did not exist 
between the in-plane and transverse equations of equilibrium. Typically, 
the accuracy of the solution of the transverse equation of equilibrium was 
strongly dependent upon the transverse deflection but only weakly influ- 
enced by the in-plane displacement. Conversely, the in-plane equation of 
equilibrium was strongly dependent upon both the transverse and i.n-plane 
deflections. Consequently, the rate of convergence of the transverse noda 
equations was much more rapid than that of the in-plane ones. 

In contrast to the treatment of the equations of equilibrium as a 
typical boundary value problem, the solution for the changes to the moment 
and force resultants between successive states of quasi-static equilibrium is 
handled as an initial value problem. As such, a fourth order Runge-kutla 
integration routine was employed to integrate the constitutive law at a 
preselected set of points across the cross section for each axial node used in 
the finite difference mesh. In this process, it was assumed that the change 
in actual stress at each of these points could be approximated as a linear 
function of time over a given, reasonably short, time interval. 

Once the “trial stresses” at the end of the time interval had been ap- 
proximated, a Newton-Cotes quadrature formula was used to numerically 
approximate the force and moment resultants. These were then employed 
to compute the deflections for this new stale of quasi-static equilibrium. 
From the deflection solution, a revised stress field could be calculated and 
compared with that which had been employed to integrate the constitutive 
law. The process was repeated if more than nominal differences existed 
between the assumed and computed changes in stress. For this, the linear 
approximating function, were adjusted based upon the computed stress dis- 
tribution. Otherwise, the results were accepted and the analysis proceeded 
on to the next time increment. 

This procedure was found to work very well after the first few time 
steps. Most of the computation effort was expended in the solution of the 
“boundary value” part of the problem and not in the iteration for material 
behavior. Principally, this was due to the fact that the rate of change in 
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actual strcsi was reasonably constant and was therefore easy to estimate 
between successive increments. This began to break down, however, when 
the behavior of the material began to resemble a plastic response. Under 
these conditions, the rate of change in actual stress would change rapidly 
even over very short time intervals. Consequently, accurately forecasting 
its rate of change was difficult. Thus, a greater number of iterations was 
needed to close the numerical loop. 

NUMERICAL RESULTS 

The problems considered are an eccentrically loaded cantilever beam- 
column and a simply supported pressure loaded shallow circular arch. In 
both cases, the behavior of the structural element is examined for constant 
loading at constant temperatures of 400, 600 and 800 °C and with sinu- 
soidal variations about the temperatures of 400 and 600 a C. Simultaneous 
variations in loading and temperature are not examined. However, the 
potential influences of time-invariant temperature gradients in the depth 
direction are examined for both constant and variable temperatures. 

The beam-column considered is 30,48 cm long having a square cross 
section of 1.27 cm depth and thickness. The eccentric load is applied 0.3048 
cm below the centroidal axis yielding an eccentricity ratio of 0.01. The 
direction of the load is assumed to remain constant. Twelve axial nodes 
are used to model the beam-column. Additionally, a five point transverse 
grid is employed at each axial node to approximate variations in stress 
and strain across the cross section. One transverse grid point is located 
at each extreme surface of the cross section and one is positioned on the 
centroidal axis. The remaining two points are spaced equidistantiy between 
these three points. It should be noted that the results from this “twelve 
axial node five grid’ 1 model compare favorably with results obtained using 
greater numbers of axial nodes and transverse grid points. 

The circular arch examined is a 8.59 deg. segment of a circle. Physi- 
cally, the arc length of the arch is 22.86 cm with an initial radius of curvature 
of 152.4 cm. For this relatively shallow arch, the rise of the centroidal axis 
is approximately 0.43 cm. Unlike the beam-column, the arch has a rect- 
angular cross section 0.51 cm in width and 0.38 cm in depth. The arch is 
divided into fourteen segments for the numerical model. Again, a five point 
transverse grid is established at the location of each axial node. Similar 
to the beam-column, the results obtained with this “fourteen segment five 
transverse grid moder compare favorably with models employing greater 
numbers of both. Illustrations of the geometry of the beam-column and 
arch models are provided in Figs. 3 and 4. 

Note that the dimensions selected for these sample problems are not 
based upon the consideration of “typical” structural elements. Instead, the 
dimensions are specifically chosen to accelerate the onset of time-dependent 
behavior. The purpose of this is to minimize the length of the initial “sta- 
ble” response period thereby reducing the overall magnitude of the compu- 
tational effort. In general, a “realistically” sized structural element would 
be appreciably stiffer and thus provide a much longer period of stable re- 
sponse. However, other than this extension of the “stable” useful life, the 
behavioral characteristics of such “realistic” structural elements would be 
highly similar to those of the example problems. 

Finally, before proceeding with a detailed discussion of the numerical 
results, some introductory remarks on the principal factors which determine 
the form of the response merit consideration. Examination of Eqn. (17a), 
for instance, reveals that the inelastic strain rate is determined by the rela- 
tive magnitudes of the actual stress and the backstress. The inelastic strain 
rate changes whenever the rate of change in actual stress varies from the 
rate of change of the backstress. In some situations, the rates of change 
of actual stress and backstress tend to equilibrate yielding a relatively con- 
stant difference. With this, the rate of inelastic straining tends to decrease 
to a nearly constant value. This behavior might be characterized as initial 
“primary creep” transitioning to “secondary” or “steady" creep . 



Figure 3 Eccentrically Loaded Cantilever Beam- Column. 



Figure 4 Axial Finite Difference Mesh for the Arch 


In other cases, the rate of change in stress increases much more 
rapidly than the backstress. Provided the magnitude of the difference 
between them is not excessive, this produces a response akin to that of 
accelerating or “tcrciary” creep. The final possibility occurs when the dif- 
ference between the actual stress and the backstress is very large. When 
this happens, the exponential nature of Eqn. (17a) creates a situation where 
inelastic strain rate tends to follow stress increment directly. Consequently, 
in the limit, the law exhibits a behavior similar to incremental plasticity. 

Thus, it should be evident that the crucial factors in the analysis 
are those which determine how rapidly the actual stress and the backsttess 
change with time. The most significant factors were found to be a geo- 
metrical effect related to the bending moment and the growth law for the 
backstress. Consequently, the overall response of the structural element 
was determined by the relative interaction between these effects. 

The geometric effect occurs in both the beam-column and the arch. 
However, it is more easily visualized with respect to the beam-column geom- 
etry and therefore is described in that context. Essentially, the maximum 
(in magnitude) stresses in the beam-column are determined by the mo- 
ment created by the end load. Transverse deflection of the beam-column, 
increasing the moment arm of the eccentric load, increases the magnitude 
of these stresses. However, tending to counterbalance this effect is the end 
rotation. Because the line of action of the load remains constant, end ro- 
tation reduces the effective moment arm created by the eccentricity of the 
load. Thus, this tends to reduce the magnitude of the end moment. 

In a prior study [25] it was found that the relative significance of these 
two influences was related to the load eccentricity and the time-dependent 
characteristics of the material. For an eccentricity of the order employed 
in this study, these two opposing effects can approximately counterbal- 
ance one another only when beam deflection is relatively minor. Beyond 
this threshold, the transverse deflection effect predominates and the end 
moment inherently increases. Therefore, the magnitudes of the maximum 
stresses increase with deflection and thus also with time. ... 

In contrast, for slowly varying temperature changes, allowing the 
0 terms of Eqn. (I7b) to be neglected, the evolution of the backstress is 
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governed by a relationship of the form 

• Oi, - (n n - n; t )6 ( 28 ) 

From Eqn. (l?e), it should evident that 6 must be non-negative. Thus, 
the magnitude of the rate of change of backstress depends upon the signs 
and magnitudes of the inelastic strain rate and the difference between the 
current and reference backstress. Because the actual stress is only indirectly 
coupled to the backstress, through the deformation and the growth laws, 
an increase in actual stress does not necessarily generate an equivalent 
increase in backstress. Thus, the quantities in Eqn. (28) need not change 
in equivalent or proportional amounts. As such, the rate of change of 
backstress may increase, decrease or remain relatively constant, thereby 
providing a wide variety of possible results. 

The Beam - Column at Constant Temperature 

The time-dependent deflection of the beam-column at constant tem- 
perature provides the simplest demonstration of these effects. Figures 5, 
6 and 7 illustrate the time-dependent end deflection of the beam-column 
under constant load at temperatures of 400, 600 and 800 • C, respectively. 
The loading is expressed in terms of the ratio of the applied load to the 
Euler load for a perfect configuration, P,, where P, =: t*.EJ/ 4L j . Note 
that the Euler load is a function of temperature due to the temperature 
dependence of the elastic modulus. Also note that the relative transverse 
deflection (i.e.: the vertical axis) represents the ratio of the time-dependent 
deflection to the initial elastic deflection. Thus, the results indicate the 
relative increase in deflection produced by inelastic straining. 

Except for the lowest loading at 400 * C, the 400 and 600 * C beam- 
column results exhibit a short initial settling period followed by a virtually 
linear increase in transverse deflection with time. This type of behavior is 
synonomous with a response of primary creep transitioning to secondary 
creep. Not unexpectedly, the higher loadings produce the greater rates of 
increase. This infers that under these conditions, the difference between 
the actual stress and the backstress must remain nearly constant with the 
greater numerical differences occurring at the highest loadings. 

This hypothesis is confirmed by in Fig. 8, which illustrates the differ- 
ence between the maximum (in magnitude) actual stress and the backstress 
for the 400 *C case.f Note that, due to the combination of bending and 
axial loading, the maximum (in magnitude) actual stress occurs in the ex- 
treme fibers adjacent to the wall on the same side of the centroidal plane 
as the applied end load. 

Except for some slight initial variations, the difference between the 
actual stress and backstress remains virtually constant. Thus, the right- 
hand side of Eqn. (17a) effectively is constant. This yields a constant 
relatively low rate of inelastic straining. Since this low rate of inelastic 
straining does not significantly alter the deflection of the beam-column, 
significant changes in the actual stresses do not occur. Concurrently, the 
low rate of inelastic straining yields a low the rate of change in backstress. 
Thus, the combination of these effects maintains an approximately constant 
difference between actual stress and backstress. 

In contrast, virtually all levels of loading at an 800 *C temperature 
produce an accelerating rate of transverse deflection. Only the lowest load 
produces a steady creep response; all higher loadings produce an accelerat- 
ing rate of deformation. At the two highest load levels, the beam-column 
deflects so rapidly it could be considered to have failed almost instanta- 
neously. Of course, in comparing results between these different cases, the 
significant differences in the time scales should be kept in mind. At the 
lower temperatures, elapsed time can be expressed in hours. At the highest 
temperature, elapsed time must be indicated in seconds. 

Why the 800 # C results are so different from those at the lower 
temperatures is directly attributable to the difference between the actual 

t Results for the difference between the maximum (in magnitude) ac- 
tual stress and the backstress at a 600 *C temperature are not included 
since they are virtually identical to those of Fig. 8 and are available else- 
where [27]. The difference tends to remain relatively constant following a 
short initial period with the greater differences at the higher loads. 
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Figures Beam-Column Deflection at 400 5 C. 



Figure 6 Beam-Column Deflection at 600 °C. 



Figure 7 Beam-Column Deflection at S00°C. 
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Figure 8. Difference between actual stress and backstress 
for the 400* C beam-column. 


stress and the backstress, illustrated in Fig. 9. At the lowest load, the initial 
relaxation is followed by a period of constant difference between the actual 
stress and backstress. Thus, the increase in actual stress from transverse 
deflection is counterbalanced by the concurrent growth in the backstress. 
Higher levels of loading alter the relative rates of growth. At higher loads, 
the rate of increase in actual stress far exceeds that of the backstress. The 
inelastic strain rate increases thus increasing the rate of deflection which, 
in turn, further increases the rate of change of actual stress. Thus, the 
process reinforces itself accelerating the approach to failure. 

Before proceeding to the arch, a few words concerning the initial 
relaxation are warranted. The “relaxation” process is a combination of 
effects. First, the inelastic deformation tends to limit the rale of increase 
in the stresses at the extreme fibers. To maintain equilibrium, load bearing 
responsibility is transfered toward the centroidal axis. Due to the relatively 
low initial magnitudes of these stresses, this normally does not produce any 
significant inelastic staining near the centroidal axis. 

The exception occurs with the behavior demonstrated in the 800 * C 
case. The rapid inelastic straining at the extreme fibers significantly in- 
creases the magnitudes of the stresses throughout the central core. Thus, 
appreciable inelastic straining occurs over the majority of the cross sec- 
tion. With this, inelastic straining along the centroidal axis also begins to 
strongly influence the overall response. In fact, the form of Tailute" which 
results might be characterised as a viscoelastic analog to a “plastic hinge. 

The Arch at Constant Temperature , 

The behavior of the pressure loaded shallow arch at constant tem- 
peratures of 400, 600 and 800 *C is generally similar to that of the beam- 
column. Figures' 10 and 11 illustrate the transverse deflection at the center 
of the arch at 600 and 800 * C. Note that the critical load for the arch is 
estimated to lie between 179 to 186 kPa (26 to 27 psi) for these tempera- 
tures. 

Again, the temperature increase from 600 to 800 # C reduces useful 
life by mote then an order of magnitude. This is apparent from the signif- 
icant difference in time scale between the figures. The differences between 
the actual stress and the backstress for these arch examples are similar to 
the those for the beam-column and thus have not been included. 


One factor .common to both temperatures is the sensitivity to load 
magnitude. 'Note that only a alight increase in pressure can substantially 
alter the character of the response. The .reason is that in the arch, the cen- 
troidal axis stress is always significant due to the curvature and boundary 
conditions. Thus, it always influences behavior. This stressing induces a 
high rate of -compressive inelastic straining along the centroidal axis “re- 
ducing” the nominal arc length of the .arch. This geometric change tends 
to accomodate additional transverse deflection through reduction of bolt, 
the “nominal” curvature and the “unloaded" arc length of the arch. Note 
also .that, unlike the beam-column, Where such an effect is localized near 
the wall, this centroidal axis .straining .occurs over most of the central sec- 
tion, of the arch. Thus, the “failure!” some tends to be distributed and not 
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Figure 9. Difference between actual stress and backstress 
for the 800" C. beam-column. 


localized. As such, the sudden rapid increase in transverse deflection which 
occurs with the beam-column is not as pronounced in the arch. 

The Influence of Time - Invariant Temperature Gradients 

The effects of a time-invariant, linear (in the depth direction) temper- 
ature gradient on the constant temperature response are considered next. 
This type of temperature gradient is associated with the steady state llow of 
heat through the depth. For simplicity, the gradient is assumed to remain 
constant with respect to the length of the element. 

For bookkeeping purposes, gradients where the temperature is great- 
est at the upper surface and lowest at the bottom surface are considered 
to be “positive.” Conversely, the case where the greatest temperature is at 
the bottom surface is denoted as a “negative” gradient. Note that a I the 
gradients were established so that the temperature of the centroidal axis 
would remain at the nominal case temperature. Thus, a +10 " C gradient 
for a 400 "C beam-column implies that the upper, centrotdal and lower 

surface temperatures are 405, 400 and 395 "C, respectively. 

This type of temperature gradient introduces two effects, irst, ue 
to the temperature dependence of the elastic modulus, the first integral of 
the Elastic modulus over a cross section does not vanish. This produces a 
weak level of bending-stretching coupling. The second effect ts that such 
gradient, cause thermal bending of the element. Positive cause 

downward bending thereby' augmenting the mechanically induced deilec 
tion. Conversely, a negative gradient reduces the load induced bending. 

It was found that neither of these exert a major influence on the con- 
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thermal gradient, the magnitude of this change is not sufficient to alter 
the occurance of the general behaviors illustrated in Figs. 8 through 11. 
These results are typical of the behavior observed for the other constant 
temperature beam-column and arch examples. 

As mentioned above, the imposition of such thermal gradients also 
creates some weak bending-stretching coupling in the governing equations. 
This was found to be an insignificant effect. Specifically, a number of test 
cases were examined where the coupling term, E u was artifically set to zero. 
There was no appreciable difference in the results with respect to results 
obtained when the coupling term was retained. Consequently, it can be 
concluded that the potential bending-stretching coupling which might be 
induced by the temperature dependence of the elastic modulus is negligible. 

Arch and Beam — Column Variable Temperature Behavior 

The final aspect of behavior examined is the response of the arch 
and beam-column under sinusoidally varying temperatures. Results are 




Figure 12. Influence of temperature difference on the deftection 
of the 400° C beam-column. 


resented for both the arch and beam-column without a temperature gra- 
ient in the depth direction. The combined influences of sinusoidal tem- 
erature variations and a time-invariant depth direction thermal grad.en 

re investigated only for the case of the beam-column. 

In the cases examined, the temperature is assumed to vary 
oida. manner about an elevated mean temperature of either 400 or 6 * C. 

tmplitudes of 50, 100 and 150 *C are employed. Also in keeping 
,se of a quasi-static analysis, 1200 and 1800 sec periods are employed 
he sinusoids. Sinusoidal variations about the 800 »C temperature were 
iot studied due to the extremely short life exhibited in the constant tern- 

P Finally, the thermal model includes a slight initial delay between 
when the load is first applied and when the sinusoidal temperature van- 
commence. With some combinations of higher 
higher loads, a step-Hke initial transient is generated in the constiti. 
law state variables. Superimposing a sinusoidal temperature at e same 
lime that this step-like increase in the stale variable, occur, forces h 
Runge-Kutta integration routine to employ exceptionally small time step 
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Figure 13. Transverse deflection as a function of temperature 
gradient for the 600® C beam-column. 

increments to provide accurate results. However, the need to employ such 
small time step increments passes almost immediately once the load in- 
duced transient change has taken place. Thus, the short delay allows the 
integration routine to stabilize before the start of the temperature oscilla- 
tions. Typically, a delay period on the order of 50 sec is employed. Based 
on a limited number of test cases, it was found that such a delay period 
has no appreciable impact on the overall results. 

Figures 14 and 15 illustrate the time dependent deflection of the 
arch and the beam-column, respectively, for 50 and 100 * C amplitude tem- 
perature oscillations about a 600 * C temperature. Note that the loading 
employed in both of these cases would nominally produce a “stable” time- 
dependent response under constant temperature conditions. Several inter- 
esting features can be discerned. First, the shape of each deflection curve is 
a distorted sinusoid. The upper peaks tend to be exaggerated whereas the 
lower peaks arc well rounded. Additionally, an underlying increasing trend 
in time-dependent deflection can be observed in both. Specifically, both 
the upper peak and lower peak deflections increase between each cycle. 



Figure 14. Center deflection for the 600* C arch for 50 and 
100* C sinusoidal temperature amplitudes. 
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Figure 15. End deflection for the 600* C beam-column for 50 
and 100* C sinusoidal temperature amplitudes. 

The distortion of the peaks from a true sinusoid is a result of the tem- 
perature dependence of the inelastic strain. During the higher temperature 
portion of the cycle, inelastic straining increases due to the temperature 
sensitivity of the constitutive law parameters. This results in a rapid in- 
crease in the deflection. Conversely, the inelastic strain rate is reduced at 
lower temperatures yielding less inelastic deformation during that part of 
the cycle. Thus, the behavior during the lower temperature portion of the 
cycle more closely approximates that of simple thermo-elastic deformation 
where deflection follows the thermal cycle exactly. 

Tables I and II summarize these increasing trends for the arch and 
beam-column cases, respectively. Note that these tables present results 
obtained for a variety of temperature differences in the depth direction. 

Table I. Hate of Change of Arch Deflection for Sinusoidal 
Temperature Variations 


Pressure 

(kPa) 

Sine 

Amplitude 

*C 

Rate of Change 

1*' to 2 nd 
Upper Peak 

(cm/sec) 

1-* to 2 nd 
Lower Peak 

165 

25 

3.18xl0" 7 

2.74xl0" 7 

165 

50 

3.84xl0" 7 

3.20xl0" 7 

165 

100 

4.98x10"* 

2.87x10"* 

152 

100 

7.87xl0" 7 

6.12xl0" 7 


As Table I illustrates, the magnitude of the increase in deflection be- 
tween the first and second upper peaks appreciably exceeds that between 
the first and second lower peaks. Although not indicated in this or the 
following table, the increase in deflection at the cycle midpoint temper- 
ature falls between the values of the extremes. It should be noted that 
this unequal expansion of the thermal cycle deflection indicates that creep 
ratchetting is occurring; 

This table also demonstrates the significant nonlinearity involved in 
this phenomena. The first doubling of the amplitude of the thermal cycle, 
from 25 to 50 * C, results in only a modest increase in the peak to peak 
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deflection*. However, a second doubling, from SO to 100 *C, increases 
the peak to peak deflection by slightly more than an order of magnitude. 
Table II illustrates that a similar effect for the beam-column, ote a 
except for the SO *C amplitude with no gradient, increasing the ampl.tudc 
of the theimal cycle increase* the change in peak to peak deflection in a 
nonlinear manner* 

Table 11 also demonstrates that the increase between consecutive 
peaks has an overall decreasing trend. The absolute increase in deflection 
decreases as the number of cycles increases. The cause of this decreasing 
trend is attributed to the growth in backstress. Generally, the maximum 
stresses do not change significantly between consecutive peaks of the ther- 
mal cycle because they are principally determined by the thermo-elastic 
deformation. The accumulated inelastic deformation does not significantly 
alter this stress distribution. However, the amount of inelastic strain is 
sufficient to create an increase in the backstress between consecutive peaks, 
one which exceeds that of the actual stress. Thus, the difference between 
the actual stress and the backstress decreases. In turn, this reduces the 
rate of inelastic straining. Of course, this effect becomes less significant as 
the backstress approaches saturation. 


Table II. 


Sine 

Ampl. 

*C 

50 

50 

50 

100 

100 

100 

150 

150 

150 


Rate of Change of Beam-Column Deflection for Sinusoidal 
Temperature Variations 


Temp. Rate Change (cm/sec) 

Diff. Upper Peak* Lower Peaks 

• C 1“ to 2 nd 2 nd to r* l ft to 2 nd 2 nd to Z' d 


0 

1.48x10'* 

1.27x10'* 

+10 

8.46x10'* 

8.26x10'* 

+20 

4.29xl0“ 7 

4.06xl0“ 7 

0 

1.07x10“* 

1.03x10'* 

+10 

2.21x10-* 

2.06x10'* 

+20 

4.72x10“* 

4.27x10'* 

0 

1.87x10“* 

1.66x10'* 

+10 

3.22x10"* 

2.91x10'* 

+20 

5.4U10'* 

4.93x10'* 


1.48x10'* 

1.27x10“* 

8.05x10“* 

7.82x10'* 

4.0U10' 7 

3.84xl0' 7 

9.83xl0' 7 

9.45xl0' 7 

2.01x10'* 

1.88x10“* 

4.24x10'* 

3.86x10'* 

1.61x10“* 

1.45x10'* 

2.77x10'* 

2.52x10“* 

4.67x10'* 

4.32x10'* 


The remaining influence illustrated by the data of Table II is the 
effect of a time-invariant depth direction temperature gradient on the in- 
elastic response. In general, the change in peak to peak deflection is signif- 
icantly greater in the presence of the time-invariant temperature gradient. 
Typically, almost an order of magnitude increase in peak to peak change in 
deflection occurs for the lowest temperature cycle amplitude. Interestingly, 
the effect of the depth direction gradient becomes less significant as the 
amplitude of the thermal cycle increases. 

The ultimate significance of this is whether or not temperature gra- 
dients in the depth direction arc truly negligible. The data of this table 
tend to indicate that the gradient may not be a negligible factor when the 
magnitude of the gradient is of the same relative magnitude as the thermal 
cycle. In such cases, the thermal gradient does appear to exert an appre- 
ciable influence on the overall response. Obviously, this also tends to imply 
that the impact of the depth direction gradient is potentially very signifi- 
cant in transient cases where the magnitude of the temperature difference 
through the depth can easily approach that of the gross transient change 
in temperature. 

CONCLUSIONS 

The elevated temperature behavior of generic types of structural el- 
ements fabricated from a typical aerospace alloy have been studied an- 
alytically using a nonlinear kinematic analysis and employing a recently 
developed nonlinear unified hereditary constitutive law to express the time- 
temperature dependence for the material. The study results demonstrate 


that, due to the specific format of the constitutive law, the behavioral re- 
sponse of the structural element is determined principally by the difference 
between the actual stress in the element and the backstress variable of the 
constitutive law. The firtt of these, the actual stress in the element, is 
basically controlled by the geometry of deformation. However, the second 
factor, the backstress, is governed by the appropriate growth rules of the 
constitutive law. 

This implies that accurate results for such an analysis can be ob- 
tained only if the form of the backstress growth law and the numerical 
values employed therein have been established to a reasonably high degree 
of certainty. Otherwise, the prediction of the conditions under which the 
response of the structural element may change from that of a “steady" form 
of creep to a rapid approach to failure could not be established with any 
degree of reliability. For entirely similar reasons, the actual stresses in the 
structural element needs to be established accurately to ensure that reliable 
results ensue. 

The specific results obtained from the sample problems indicate that 
a. constitutive law of the type proposed by Walker has the capability to 
model various forms of creep behavior. Under various constant load and 
temperature conditions, it has been shown that the predicted response may 
vary from simple primary creep-secondary creep to almost instantaneous 
terciary creep. Additionally, under varying temperature conditions, the 
Walker law has demonstrated the capability to predict the development of 
the type of inelastic strain biasing which can produce creep ratchetting. 

Assuming the Walker law and associated constants provide an accu- 
rate representation of the elevated temperature response of a material such 
as Hastelloy X, the sample problem results indicate that the existence of a 
depth direction temperature gradient in n thin section may not always be 
negligible. Specifically, under constant temperature conditions, the results 
indicate that such temperature gradients do not appear to exert any sig- 
nificant influence on behavioral response. However, such gradients do have 
a significant impact when the. overall tempernture of the element is not 
constant. The existence of this type of gradient can appreciably accelerate 
the overall rate of creep. It is indicated that this specific effect is most 
pronounced when the magnitude of the temperature difference through the 
depth of the element is on the order of the amplitude of the gross variation 
in temperature of the element. Consequently, in addition to when such con- 
ditions may exist in a “quasi-steady state" situation, this implies that such 
gradients may exert considerable influence in transient thermal problems. 
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APPENDIX A CONSTITUTIVE LAW 

In its most general integral form, Walker’s functional theory has the 
format: 


<*,(<).= + Sa J‘ (A[©(£)] + |m(©(01) - 3<*[©«)1 J) 

(A ~ Id) 

fM«) = n«(‘) + *i(©(<)M‘) + «*(©(»)] ^ di, 

(4-16) 

K(l) = K, [©(<)]- K,[©f$e— ’ r[*(*)|w(«), (.4 - lc) 




*i 

B£ 


d<Tij 

Bi 


-<y O (6(O)(3A[0(«)l + 2/l[©({.)]) J) M, 




nf, = -0‘ 


Cq(t)c ti («) 

9 


). 


(.4 - Id) 


(A - U) 


run - /" 3/*l©(0] , 9W .i-i/.[e«)) 
WJ- y. K(f) ( 6i> 


G(i) - J'{ (n,(©(C)l + „ 4 (©(C)]e-< G «)l w «>) ^ 


(.4 - »/) 


(A - 1 g) 


and 



y 3 Bi Bi ^ 


(.4 - Ih) 


Square bracket terms, of the form /*[©.(<)] , are used to denote the depen- 
dence of the material constants A, /*, fl*, n, m, n t , nj, nj, n*, n*, n*, n 7 , 
Ki, and Kj on the temperature, O. However, for clarity, the indication of 
the explicit dependence on temperature will be suppressed in the following 
development. 

The differential format for Walker’s functional theory is obtained 
through differentiating the above relations * : th respect to time. Employing 
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Leibnis’s rule, differential formats for Eqns. (A-ld), (A-lf), (A-lg) and 
(A-lh) are readily obtained. A somewhat simplified differentia! format 
for Eqn. (A-lb) may be obtained in the following manner. Differentiating 
Eqn. (A-lb) with respect to time and using Leibnis's rule, after rearranging, 
yields 

A„ - Af, +(., + |&* - .,61 

(>1-2) 

The integral which appears in the above is identical to the integral of 
Eqn. (A-lb). Consequently, solving Eqn. (A-lb) for the integral and substi- 
tuting into Eqn. (A-2) yields the differential format given by Eqn. (A-6b). 

With some manipulation, the differential form of Eqn. (A-la) may 
be stated in a “power law” type format. First, Eqn. (A-la) is differentiated 
with respect to time and then simplified using the deviatoric stress tensor 
to yield 


i(f« y - m = jf * 

(A - 3) 

Note that to establish the above it is necessary to show that flu = flj k = 0. 
These, however, follow directly from Eqns. (A-lb) and (A-le) provided the 
inelastic portion of the deformation may be treated as, at least, approxi- 


mately incompressible (i.e., c*i « 0). 

Therefore, substituting Eqn. (A-3) into Eqn. (A-2) and using the 
differential form of Eqn. (A-ld) yields: 


ij - — Q ( 2 S * J ~ ^ 

The Q term can be replaced with the aid of the differential forms of 
Eqns. (A-lf) and (A-lh). After some algebraic manipulation and minor 
rearranging this yields: 

w = (. 4 - 5 ) 

Consequently, substituting Eqn. (A-5) hack into the differential form of 
Eqn. (A-lf) and using Eqn. (A-4) yields, after some algebra, the differential 
power law of Eqn. (A-la). Thus, the complete set of differential relations 
for Walker's functional theory are: 


iii = { 


- n*i)(|su - n *i) \ 


v/hF 


(fsjj fl»j) 


n) (| S »nn *“ Amn) 

{A - 6a) 


n<, = n*; + (m + _ ( n 'i ~ n *j - n,Ci >) (® - fie ®) 1 

(A - 66) 

K = K, - K,e-"' w , (A -6c) 


2/*c ij = SijXihk + 2/ic ij — arij — £»j( 3A + 2ft)a©, (A — fid) 


p* t ^Q# 3fl* 1 

fl?. = Trr + ; „ \i - 2ctsCkjC Ff c pf ], 

(A - fie) 


and 

W = (A-6 9 ) 

Note that a differential form for the function Q is not required since it does 
not appear in any of the other expressions. 

Numerical constants for the law were established [22-24] from uniax- 
ial bar specimens tested under fully reversed, strain-rate controlled cyclic 
stress-str; in tests. These tests were conducted for a variety of temperatures 
and strain rates. The testing conditions were sufficiently rigorous to cause 
plastic deformation during the loading cycle. The creep and relaxation 
‘ properties were deduced from observation of the behavior of the samples 
when the cyclic loading was “held” at various points on the “steady-state 
hysteresis loop. Appendix Figs. A-l through A-7 illustrate the tempera- 
ture dependent parametric values established under that work. Note that 
the material constants K 2| n t , n 4l n», and n are sero over the enitre tem- 
perature range considered. This results in additional simplifications to the 
constitutive law. 

APPENDIX B FINITE DIFFERENCE EQUATIONS 

Solutions for the governing nonlinear differential equations are ob- 
tained using a Newton type method for nonlinear differential equations. It 
is assumed that trial solutions for the transverse and axial deflection func- 
tions, denoted as w and u, respectively, exist. From these, trial values for 
centroidal axis strain and cross section rotation, e 0 *nd respectively, are 
determined. Then, using the approximations e Q ~ €„ + At. and <t> = <j>+&4> 
to substitute into Eqns. (16) yields 


&±<f> 0 a AM 


I(E 0 A* # +B l ^)-^F = P*-JtN + 


0 a M 

0, a * 


(5 -la) 


(Eo + iEO^-M^ + iE*) 


d'A<t> 

d. 7 


an -*M 

8s ~ * ds ' 


(B - 16 ) 


The terms N and M are evaluated by substituting and 4> into Eqns. (ll). 
Numerical method? are used to approximate the derivative terms on the 
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Figure A-l. Young’s modulus as a function of temperature 
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Figure A-2. K\ as a function of temperature. 



right-hand side of Eqns (B-l). The second derivative of the moment cor- 
rection, which appears on the left-hand side of £qn. (B-la), also is evalu- 
ated numerically. This avoids the numerical problems which Can result in 
approximating fourth derivatives directly. 

Assuming w = w + Aw and a similar expansion for u, upon sub- 
stitution into Eqn. (7) and neglecting terms of order (Aw)* and higher, 
yields 

A<. = A(^ + «Aw) + B(^-*Au) (B-2) 

where the factors A and B are defined as: 



Figure A-4 n : as a Function of Temperature 



Figure A- 5 n 3 as a Function of Temperature 



Figure A-6 rt 6 as a Function of Temperature 
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Figuie A-7. Reference backstress, 0*, as a function of temperature. 


A =s 1 + — + kw and B = — - «u. (B - 3a, 6) 

8s os 

Differentiating Eqn. (B-2) with respect to s yields, after rearranging, 

*>• S, 0. (B_4) 

+B £1^ + (D + kA)£^ + kCAw 
where C and D are given by 


C = 


d 7 n 8w 

— 7 4- K-r- 
d* 1 Os 


and 


_ _ 

“ 8s 1 *8s' 


(B 5a, 6) 


Developing an expression for A<£ is accomplished by using Eqn. (8). Since 
it is assumed that <f> — $ 4* A^, substitution into Eqn. (8) yields 


—flw 

8s 


4- 


— 5Aw 

ds 


4- kAu 


**+ + ^ * + A * C05 * * \/l + 2<. + 2a 7. + ^l + W. + JA*. 

/ o _ c 


(B-6) 


where the approximations sin A </» fts A<£ and cos A ^ w 1 have been used. 
For small <. it is reasonable to assume that Ac. < e., thus allowing the Ac. 
term to be neglected. Thus, since the first term on the left-hand side, 
according to Eqn. (8), is equal to the first term of the right-hand side, this 
implies that 


— + kAu 

A</> COS <fi % ** ’ -■:=■■ . (# “ 7 ) 

y/l 4- 2c. 

Therefore, combining Eqns. (8) and (B-7) yields 

A * w x ( "^r +, ‘ Au)- (»-«) 

Derivatives of A^ are obtained from Eqn. (B-8). In this, it is assumed 
that the term 1 /A may be treated as an approximate constant and thus 
need not btf differentiated. This provides significant simplifications without 
engendering any substantial inaccuracies since the correction terms become 


negligible as the exact solution is approached. Consequently, differentiating 
Eqn. (B-8) with respect to s yields: 


and 


8A<I> 1. fl 5 Aw 0Au v 


(B - 9a) 


fi’A* 1, 0* Aw 0*Au, 

~B?~ " A * 57“ + 


(B - 95) 


Central difference formulae are used to approximate the derivative 
terms. However, prior to substituting into Eqns. (B-l) it is useful to estab- 
lish the following definitions: 


bi = kjEi + N< 4* — \ , (B — 10a, 6) 
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These two equations provide the finite-difference approximation for the de- 
flection of the internal nodes of the beam, that is for the subscript » over 
the range 3 < i < n - 1, (see Fig. 4). Therefore, for an arch of n + 1 
total nodes, they provide a system of 2n - 6 equations in 2n + 2 unknowns. 
Thus, eight additional equations are required to provide a unique solution 
for the problem. Six of these equations are obtained by using the boundary 
condition. These provide two sets of three equations, one set applicable to 
node 1 and the other set for node n + 1. 

As evidenced above, the boundary conditions do not provide a suffi- 
cient number of relations to enable unique solution. This, in part, results 
from the inherent coupling generated by the assumption of Eulet-Bernoulli 
bending when both axial and transverse deflections are possible. Note that 
the deflection functions are not only related to one another, but they also 
appear in the same functional format in both expressions. The derivative 
of one is added to (subtracted from) the other multiplied by the curvature. 
This, in turn, indicates that a form of implicit coupling exists between the 
two generalised displacements, centroidal axis strain and rotation. Ideally, 
the generalised displacements should be completely independent. 

A number of approaches can be employed to generate the two addi- 
tional equations needed for unique solution. The one used in this study is 
to require the “centroidal axis" strain in the wall to vanish. Mathemati- 
cally. this is equivalent to appending, to the existing system of equations, 
the two additional equations: 

Ae., = -?.i and Ac„ = -?«. (B-14a,i) 

where and t„ represent the strain in the left and right extensions, respec- 
tively. No condition is established for the “rotations* which might occur in 
the wall sections. This is not felt to be significant since the nodal mesh, at 
least for the equilibrium problem, is a centroidal mesh. 

Using the boundary conditions and the two “wall sttain” relations it 
is possible to eliminate the transverse and axial displacement components 
for nodes 1, 2, n, and n + 1 from the general finite-difference equations. 
This leaves a system of exactly 2n - 6 equations in the same number of 
unknowns. If the coefficients of these equations ate arranged in matrix 
form, the result is a banded matrix with a bandwidth of six. 


56 



APPENDIX B 
Novel Approach Papers 


57 



Thermodynamically Consistent Constitutive Equations 
for Nonisothermal Large-Strain, Elastoplastic, Creep Behavior 

Richard Riff,* Robert L. Carlson, t and George J. SimitsesJ 
Georgia Institute of Technology, Atlanta, Georgia 


The paper is concerned with the development of constitutive relations for large nonisothermal elastic- 
viscoplastic deformations for metals. The kinematics of elastic-plastic deformation, valid for finite strains and 
rotations* is presented. The resulting elastic-plastic uncoupled equations for the deformation rate combined with 
use of the incremental elasticity law permits a precise and purely deductive development of elastic-viscoplastic 
theory. It is shown that a phenomenological thermodynamic theory in which the elastic deformation and the 
temperature are state variables, including few internal variables, can be utilized to construct elastic-viscoplastic 
constitutive equations appropriate for metals. The limiting case of inviscid plasticity is examined. 
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= element of area 
= material line element 
= deformation rate 
= strain rate 
= yield function 
» deformation gradients 
= base vectors 
= metric tensors 

^absolute determinant of the deformation 
gradient 

= internal variables 
- normal to the surface 
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= specific applied heat 
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= temperature 
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= time 
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= density 

= Cauchy stress tensor 
= Kirchhoff stress tensor 
= specific free energy 

= Jauman stress rate 
= time derivative 
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Introduction 

T HE prediction of inelastic behavior of metallic materials 
at elevated temperature has increased in importance in 
recent years. Many important engineering applications in- 
volve the use of metals subjected to cyclic thermomcchanieal 
loads, e.g., hot section components of turbine engines, nu- 
clear reactor components, etc. These materials exhibit 
substantial complexity in their thermomechanical constitu- 
tion. In fact, so complex is their material response that it 
could be argued that without useful a priori information, ex- 
perimental characterization is futile. It is, therefore, impor- 
tant to be able to model accurately the nonelastic behavior 
of metals under cyclic mechanical and thermal loading at 
temperature levels for which creep and recovery introduce 
significant response phenomena. 

Under this kind of severe loading conditions, the real 
world of structural behavior is highly nonlinear due to the 
combined action of geometrical and physical nonlinearities. 
On one side, finite deformation (in a stressed structure) in- 
troduce nonlinear geometric effects. On the other side, 
physical nonlinearities arise even in small strain regimes, 
whereby inelastic phenomena play a particularly important 
role. From a theoretical standpoint nonlinear constitutive 
equations should be applied only in connection with 
nonlinear deformation measures. However, in engineering 
practice, the two sources of nonlinearities are separated for 
practical reasons, yielding at one end of the spectrum large 
displacement and large rotation problems and on the other 
end inelastic analysis in the presence of small strain. 

Constitutive models for small strain in engineering litera- 
ture may generally be grouped into three categories: 
classical plasticity, nonlinear viscoelasticity, and theories based 
on microstructural phenomena. Each group can be further 
separated into “unified” and “uncoupled” theories, where 
the two differ in their approach to the treatment of rate- 
independent and rate-dependent inelastic deformation. The 
uncoupled theories decompose the inelastic strain rate into a 
time-independent plastic strain rate and a time-dependent 
creep rate with independent constitutive relations describing 
plastic and creep behavior. Such uncoupling of the strain 
components provides for simpler theories to be developed, 
but precludes any creep/plasticity interaction. Recognizing 
that cyclic plasticity, creep, and recovery are not independent 
phenomena but rather are very interdependent, a number of 
“unified” models for inherently time-dependent nonelastic 
deformation have been developed recently. 
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Classical incremental plasticity theories are macro- 
phenomenological because they base the derivation of state 
variables purely on experimental results without direct 
reference to the microstructure of the material. Most in- 
cremental plasticity theories have four major components: 1) 
a stress-elastic strain relation, 2) a yield function describing 
the onset of plastic deformation, 3) a hardening rule that 
prescribes the strain hardening of the material and the 
modification of the yield surface during plastic flow, and 4) 
a flow rule that defines the components of strain that is 
plastic or nonrecoverable. Research in this area is volumi- 
nous. For example, Zienkiewicz and Cormeau 1 developed a 
rate-dependent unified theory that allows for nonassociative 
plasticity and strain softening, but does not model the 
Bauschinger effect or temperature dependence. Extensions of 
classical plasticity to model both rate and temperature effects 
were presented recently by Allen and Haisier, 2 Haisler and 
Cronenworth, 3 and Yamada and Sakurai. 4 

In the nonlinear viscoelastic approach, the constitutive 
relation is expressed as a single integral or convoluted form. 
This type of constitutive model employs the thermodynamic 
laws along with physical constraints to complete the for- 
mulation. A detailed review of several existing theories is 
presented by Walker. 5 Walker’s theory is based on a unified 
viscoplastic integral developed by modifying the constitutive 
relations for a linear three-parameter viscoelastic solid. The 
theory contains clearly defined material parameters, a rate- 
dependent equilibrium stress, and a proposed multiaxial 
model. An important shortcoming of Walker’s theory is its 
failure to model transient temperature conditions. Many 
other nonlinear viscoelastic theories have been proposed, in- 
cluding those by Cernocky and KrempI, 6 Valanis, 7 and 
Chabache. 8 

The microphenomenologicai theories attempt to represent 
the response of polycrystalline materials in terms of various 
micromechanisms of deformation and failure. Various 
dislocation theories have been developed to predict plastic 
deformation in terms of dislocation interaction, slip, glide, 
density, etc. Most of the material models developed to date 
depend primarily on the number of state variables used and 
their growth or evolutionary laws. Many of the recent 
“unified” microphenomenologicai theories have been dis- 
cussed and evaluated by Walker 9 and Chan et al. 10 

One example of a microphysically based constitutive law is 
an elastic-viscoplastic theory based on two internal state 
variables as proposed by Bodner et al. 1 * These ' authors 
demonstrate the ability of the constitutive equations to repre- 
sent the principal features of cyclic loading behavior, in- 
cluding softening upon stress reversal, cyclic hardening or 
softening, cyclic saturation, cyclic relaxation, and cyclic 
creep. One limitation of the formulation though is that the 
computed stress-strain curves are independent of the strain 
amplitude and therefore too “flat” or “square.” 

Miller 12 has reported research on the modeling of cyclic 
plasticity with “unified” constitutive equations. He also 
recognizes the shortcomings of many theories in predicting 
hysteresis loops that are oversquare in comparison to observed 
experimental behavior. Improvement is accomplished by 
making the kinematic work-hardening coefficient depend on 
the back stress and the sign of the noneiastic strain term. 
Theories that are similar in format to Miller’s have been pro- 
posed by Krieg et al. 13 and Hart. 14 The models use two inter- 
nal state variables to reflect current microstructure state and 
are based upon models for dislocation processes in pure 
metals. All of these constitutive theories were formulated 
without the use of a yield criterion. Since these models do 
not contain a completely elastic regime, the function that 
describes the inelastic strain rate should be such that the in- 
elastic strain rate is very small for low stress levels. Theories 
with a yield function and a full elastic regime have been 
developed for the case of isotropic hardening by Robinson 15 
and Lee and Zavrel 16 for both isotropic and directional 
hardening. 


As previously noted, the quantities utilized in the small 
strain theory of viscoplasticity (stress, strain, stress rate, and 
strain rate) are defined only within the assumption of “small 
strains.’* Yet the precise definition of what constitutes 
“small strain” is always left unstated. Whether or not the 
stresses for a given case are “small” cannot be determined a 
priori by geometric considerations. In general, one cannot 
know in advance whether, for a given loading of a material, 
the “small-strain” assumption (always left undefined) will 
hold or not. The question of whether the small-strain ap- 
proximations are valid is always avoided in the “small- 
strain” literature. Furthermore, as Hill 17 points out, the really 
typical plastic problems involve changes in geometry that 
cannot be disregarded. In many cases, for example, it is suf- 
ficient to take into account finite plastic strains and small 
elastic strains or vice versa. From the theoretical viewpoint, 
it is desirable in all cases to have a theory that intrinsically 
allows for both the elastic and plastic strains to be large. 
Such a theory, of course, must reduce to the earlier mentioned 
special cases, as limiting cases. Furthermore, such theories 
provide a check for those obtained by generalizing small- 
strain theories. 

The mathematical theories of deformation and flow of 
matter deal essentially with the gross properties of a 
medium. Heat and mechanical work are considered as addi- 
tional causes for a change of the state of the medium. The 
resulting phenomena in any particular material are not 
unrelated. Therefore, a thermodynamical treatment of the 
foundation of the theory of flow and deformation is ap- 
propriate and, indeed, the obvious approach. Two very dif- 
ferent main approaches to a thermodynamic theory of a con- 
tinuum can be identified. These differ from each other in the 
fundamental postulates upon which the theories are based. 
An essential controversy (a good survey of this controversy is 
given in Ref. 18) can be traced through the whole discussion 
of the thermodynamic aspects of continuum mechanics. 
None of these approaches is concerned with the atomic struc- 
ture of the material. Therefore, they represent purely phenom- 
enological approximations. Both theories are characterized 
by the same fundamental requirement that the results should 
be obtained without having recourse to statical or kinetic 
methods. 

Within each of these approaches, there are two distinct 
methods of describing history and dissipative effects: the 
functional theory 19 in which all dependent variables are 
assumed to depend on the entire history of the independent 
variables and the internal variable approach 20 wherein 
history dependence is postulated to appear implicitly in a set 
of internal variables. For experimental as well as analytical 
reasons, 21,22 the use of internal variables in modeling in- 
elastic solids is gaining widespread usage in current research. 
The main differences among the various modern theories lie 
in the choice of these internal variables. 

Therefore, the predictive value of an elastic-viscoplastic 
material model for nonisothermal, large-deformation anal- 
yses depends on three basic elements: 1) the nonlinear 
kinematic description of the elastic-plastic deformation, 2) 
thermodynamic considerations, and 3) the choice of external 
and internal thermodynamic variables. The objective of this 
paper is to examine each of these elements, illustrate their in- 
teraction, and extend these considerations to model the 
large, nonisothermal, elastic-viscoplastic deformation 
behavior of metals. 

Moreover, the paper deals with the phenomenological 
theory of elastic-viscoplastic bodies. The process inside the 
lattice and at the border of the crystal grains is taken as the 
physical background, without considering its connection to 
the macroscopic behavior of the material at the present. 

Kinematic and Fundamental Considerations 

Consider body of volume V that occupies a finite region 
of Euclidean space. When subjected to prescribed body 
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forces, surface tractions, surface temperature, and surface 
velocities, the body undergoes motion characterized by 
x i -x i (X°,t). The material particles of the body are iden- 
tified by coordinates X a , which are referred to as material 
coordinates. The relation of the material particles to the 
material coordinates X a does not change in time- The places 
in space that the particles occupy during the motion are iden- 
tified by the coordinates x i . Functions x‘ describe the motion 
of the particles X a through space. The place occupied by the 
body at * = 0 is taken as the initial configuration, in this con- 
figuration the body is assumed to be strain-free, but not 
necessarily stress-free. 

A third coordinate system is defined by the material coor- 
dinates as they deform with the body. This system will be 
denoted by X A * which are referred to as converted coor- 
dinates. The current configuration of the body with spatial 
coordinates x i and convected coordinates X A and the initial 
configuration of the body with material coordinates X a will 
be employed in what follows. For the spatial coordinates x* % 
the covariant base vectors g rt the contravariant base vector 
g r t the metric g„, and its dual g n are used. Similarly, for the 
convected coordinates X A t the covariant base vectors G A , 
the contravariant base vectors G A , the metric tensor G AB , 
and its dual G AB are used. With regard to the initial con- 
figuration, the covariant base vectors Ga, the contravariant 
base vectors G*, the metric tensor G a $, and it dual are 

used for the material coordinates X a . 

For a second-order tensor A with components A n in the 
spatial coordinates and components A AB in the convected 
coordinates, the following is true: 

A -A n g r g s =A ab G a G b (1) 

The two sets of components are related to each other 
through 

A**~* a x S ' B A ab (2) 

where X r A denotes the partial derivative dx r lX B ,t)/dX A . 

For the motion, characterized by x r (X A J) = \ r ( 1 ) > we 
have 

Ga-^.aSt Gab + ^.aX^bSts ( 3 ) 

From Eq. (3), it is seen that G AB = 0, where the dot denotes 
time material derivative. The tensor transformation equa- 
tions (1) and (2) will be used extensively in what follows. 

A material line element ds = dX a G a in the initial con- 
figuration when subjected to motion x r {X*J) is deformed 
into ds-dx f g f in the current configuration. The line element 
dxf is related to the line element dX a through the deforma- 
tion gradient F r a by dxf - F r Q dX a where 

(4) 

The mapping defined by the deformation gradient F= 
F r n g r G ' allows one to shift quantities from the current con- 
figuration to equivalent, but alternate, quantities in the in- 
itial configuration. For example, the right Caushy-Green ten- 
sor C = C alS G a G‘ J and the Green-St. Venant strain tensor 
E = E tti G a G* t in the initial configuration are 

ds = dS = g rs dx r dx s - G a jdX«dX 3 

-(grsF^-G^dX'dX* 

= ( C ni} -G Qti )dX«dX* = 2 E a3 dX°dX* (5) 

The components of the deformation gradient, which relate 
a deformed line element dX A in the convected coordinates to 
the undeformed line element dX" in the iffliriaf configuration. 


are given by j A > 

( 6 ) 

Equation (6) places in a single expression the easily confused 
but distinct ideas of the transformation of tensor com- 
ponents under a change of coordinates and a shift between 
the current configuration and the initial configuration as a 
setting for the governing equations. Truesdell and Toupin 21 
and Truesdell and Noll 24 emphasizes the current configura- 
tion with the spatial coordinates and an initial configuration 
with material coordinates. As a result, the deformation gra- 
dient plays a prominent role in their work. Only in isolated 
spots do they mention convected coordinates and, then, as 
indirectly as possible. On the other hand. Green and 
Adkins 23 and Sedov 26 rely heavily on convected coordinates. 
Our intention here is only to tie the two together for the pur- 
pose of discussing elementary assumptions. Recently, Men- 
delssohn and Baruch 27 review this same point as well as ad- 
ditional material relevant to sound numerical formulation of 
finite deformation problems. 

The velocity v = tfg r of a particle X" is defined by 

(7) 

From the spatial gradient of the velocity, the deformation 
rate 

d = d rs g r f = d AB G A G H (8) 

is defined as 

drs~ X MK.s + Kr) ( 9 ) 

The spin 

W=W n g r tf=W m G A G H (10) 

is defined as 

^rs^ l MK.,-K.r) (ID 

In the initial configuration, the Green-St. Venant strain rate 
is the shifted deformation rate, 

E<*=nF'id n =ttttd AB (12) 

Basic to most of the postulated models of large elastic- 
plastic deformation behavior is the additive decomposition 
of d rs and E n3 into elastic and plastic parts, 2S 

£ P £ P • 

d n ** d n + d n , E uli — E, vi + £"$ (13) 

The validity of this additive decomposition in the case of 
finite elastic-plastic strains has been questioned by Lee and 
his associates. 29 ' 32 Lee’s 29 approach is based on the total 
purely elastic unloading from the current state to an in- 
termediate unstressed plasticity deformed configuration, 
without any reverse or other kind of plastic flow. The major 
point in his theory is to decouple the total elastically induced 
distortion and measure it from a relaxed unstressed state, 
which is only plastically deformed from the initial to the in- 
termediate configuration. Accordingly, the deformation gra- 
dient F is decomposed in the form 

EP 

F=F F (14) 

P 

where F transforms a line element from the initial configura- 

% E 

lion to the intermediate configuration and F from the latter 
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to the current configuration. The intermediate configuration 

P # 

is chosen in such a way so that F is unaffected by the 
presence of rigid-body motion. The deformation rate tensors 
d sr and d„ are then defined. After some manipulations, Lee 
shows the following relation: 

E E P E I E P I 

d rs =<T + F rk d kt F~ i | + F rk \Y r J h ' | (15) 

where the subscripts s denote the symmetric parts. Gen- 
eralization of Lee’s theory for anisotropic elasticity was 
given by Mandel. 33 

Lee’s theory is based on the assumption that the elastic 
law does not change with the history of deformation and, 
hence, a total elastic unloading can take place. However, it 
has been shown 34 that after a fair amount of plastic flow has 
taken place, reverse plastic deformation will result soon 
upon unloading, even for small strains. Therefore, a total 
elastic unloading cannot have any physical significance. In 
view of this, the theory of Lee appears as a special case of 
the theory of Green and Naghdi. 35 Although not as general 
as the theory of Green and Naghdi, Lee’s theory has the ad- 
vantage of being more easily fitted with the physical property 
of invariance of elasticity with respect to plastic deforma- 
tion. In particular, Mandel 33 has pointed out that the Green- 
Naghdi theory is not convenient if one wants to include 
anisotropic elasticity effects. AH this can be avoided by the 
use of the convected coordinates, as proposed by Sedov 26 
and Lehmann. 36 The formulation presented herein will 
follow the work of Lehmann, 

All quantities from here on will be related to the metric of 
the coordinate system X A in the deforced state. Hence, 


A^G'"*G hC 

= (16) 


and the deformation rate is 


d A — ViG AH G - - ViG ca G*" 

= < 17 > 

The deformation gradient may be split into its elastic and its 
plastic components in the following manner: 

P E 

F\ b A 

(18) 

E P 

1 )r 


The use of capital greek subscripts and superscripts (G jr ) 
denotes parameters belonging to a fictitious intermediate 
state, which is in general incompatible. The circumstance of 
the noncontinuous configuration in the unstressed state has 
been observed by Sedov, 26 who points out that convected 
coordinates, as used herein, become non-Euclidean in this 
configuration. 

This multiplicative splitting of the metric change in the 
convected coordinates leads to an additive splitting of the 


deformation rate according to 

E E PE 

d* = sy m 14 I (/- 1 )(! (/) r f I + sym '/: | (/-')?(/)* r /. r I 


E E EP 

= sym'/ 2 | -sym 1 /:) (/ ‘)r (/” 1 ) r s/c I (19) 


E P 

d£ + d£ 

In the current configuration of the body Y t consider an 
element of area da on the surface of 5 with an outward nor- 
mal n = n r g r - n A G A . If the force d P=dP r g r = dP A G A is act- 
ing on this element, the fraction vector is / = dP/da, The 
Caushy stress, 

a = o^grSs -o AB G a G b (20) 

defined, such that f = <w i, which in component form (in 
terms of spatial coordinates) is In the convected 

coordinates, it is t A =o AB n B . 

It is convenient to work with the Kirchhoff stress tensor r 
in the current configuration, obtained from the Cauchy 
stress by scaling 

tb = — °bJ°b ( 21 ) 

P 

where p denotes the current mass density, p () the mass density 
in the initial state, and J the absolute determinant of the 
deformation gradient at the current configuration. 

The time derivative of a tensor such as stress, which is 
associated with the current configuration, admits infinitely 
many definitions, depending upon the coordinate system 
employed in the computation of such time derivatives. For a 
correct large-strain, large-rotation elastic-plastic model, the 
notion of “invariant stress fluxes” and “objectivity” must 
be introduced. A good treatment requires more space than is 
available here. 23 The corotational stress rate, also referred to 
as Jauman stress rate, will suffice for the purpose of the 
present discussion. Hence, in convected coordinates, 

v 

Og =5 . 4 * d^Cg — dgO’c 

^ + dcTg — dgTf (22) 

From Eq. (21), the following relations between the various 
rates of Kirchoff stress and Cauchy stress are obtained: 

*i+Jd c c oi 

P 

^ B = ^i + Jd c c <ri (23) 

P 

If a rate constitutive law is postulated between a and d in 
finite inelasticity theories, then a potential does not exist, 
which is necessary in the variational or thermodynamics- 
based formulation of the problem. The basic difficulty lies 
with the d f t . term. This is remedied by postulating a con- 
stitutive law between r and d. 

The Elastic Deformations 

The present study is concerned with the structure of the 
constitutive relation of an elastic-viscoplastic (elastic-plastic) 
medium. The term elastic-viscoplastic means that the viscos- 
ity does not intervene in the elastic domain whose boundary, 
in particular, is well defined at every stage of the deforma- 
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tion. For simplicity, we assume further that the thermoelastic 
behavior of the body is isotropic and unaffected by inelastic 
deformation in the sense that the material constants 
characterizing the thermoelastic behavior are independent of 
inelastic deformation. Thus, we can obtain a unique relation 

£ 

between the elastic deformations represented by /£, the Kirch- 
hoff stresses r£, and the temperature 7\ 35 - 36 

£ £ £ £ 
fc = /c( r c» T)\ r^-T^trfcr); T-T(rtJt) (24) 

This function may be transformed into an incremental rela- 
tion by differentiation with respect to time. This leads to a 
general expression of the form 

£ £ V 

dt = d£\rt,r£,T,t,G AC ,d£) (25) 

From Eq. (24) we see that the total deformation rate enters 
the incremental form of the thermoelastic stress-strain rela- 
tions. Therefore, the thermoelastic deformation is not in- 
dependent of the inelastic deformation occurring at the same 
time. This follows from the fact that in the integrated form 
of the thermoelastic stress-strain relations Eq. (24), the 
stresses and the strains are referred to the deformed state of 
the body. 

In view of the present discussion and the discussion in the 
previous section, the hyperelastic behavior described by Eqs. 
(24) and (25) will be replaced by a hypoelastic law. The 
hypoelastic law is a path-dependent material law, since it 
cannot be expressed in terms of an initial and a final state; it 
depends on the path connecting these states. Otherwise, if we 
did not make such a change, it would be necessary to retain 
the finite deformation measure in the constitutive law. For 
small elastic strains, there is practically no difference be- 
tween hypoelastic and hyperelastic laws, as shown, for exam- 
ple, by Lehmann. 36 

The above could be illustrated by the following example. 
From the frequently used elastic stress-strain relation, 

* = ‘ /J 1 ' : *} = 1g y - 777 T?5 ") 


temperature. The inequality given by Eq. (29b) is the loading 
condition. The actual form of the yield condition for a given 
material is determined by a set of so-called internal 
parameters, which are scalars and/or tensors of even order. 
The current values of the internal parameters depend on the 
initial state of the material and the history of the ther- 
momechanical process. 

Thermodynamic Processes 

In the treatment of elastic-plastic or elastic-viscoplastic 
deformations, we have to distinguish between the description 
as a thermomechanical process and the corresponding one by 
means of thermodynamic state equations. It is sometimes 
assumed that, in the case of processes which proceed through 
nonequilibrium states, it is fundamentally necessary to start 
with a description of the process. 19 - 24 ’ 37 Alternatively, it has 
been proposed that one might assume local equilibrium for 
the elements of a body and therefore describe the state of the 
elements, in general, by state equations. 36 " 40 The conse- 
quences of adopting these two approaches become par- 
ticularly clear when considering the influence of entropy. In 
the description of the process, entropy is a derived quantity 
and in principle we can proceed without introducing it. In 
the description by state equations, it is, on the contrary, a 
necessary state value that, at least in principle, can be im- 
mediately determined. When restricting ourselves to homo- 
geneous, quasistatical thermomechanical processes, the 
description by state equations can be reviewed as equivalent 
to that by processes. 37 - 41 The controversial issues will, thus, 
not be discussed further. 

Restricting ourselves to elementary processes, we need not 
analyze whether the applied heat arises from heat conduction 
or from heat sources. For the same reason, it is not 
necessary, in our case, to introduce the temperature gradient 
in addition to the temperature or the body forces in addition 
to the stresses. 

The first law states, under our simplifying assumptions, 
that the rate of the specific internal energy u is the sum of 
the rates of the specific mechanical work w and the specific 
applied heat q , 

u = w + q (31) 


+ a T 0 )6'c 


(26) The rate of mechanical work is given by 


we get 

£ i r r£ 1 v £ -) £ 

</c= 2^ + «77c (27) 

which may be replaced by 

E 1 f 7 v v 

^ = 2cN" 7T7 rf^J+o7^ + «7^ (28) 

We assume that inelastic deformation occurs if and only if 

F(r A c ,T,k a A c . A A C B D ,...) =0 (29a) 

BF v dF . 

_* + _7>0 <»b> 

or, for elastic-plastic material, 

F(rt,T.k a£,...,A£g,...)>0 (30) 

for an elastic-viscoplastic material. The function £ represents 
the yield condition that bounds the domain of pure thermo- 
elastic behavior in the ten-dimensional space of stress and 


Po 


(32) 


and may be split into an elastic and an inelastic part accord- 
ing to Eq. (19), 


i £ t PEI 

W= — r£d c A + — r$(f A = W+ W (33) 

Po Po 


D 

The rate of inelastic work must also be split into a part W , 

which is dissipated at once, and into another part W % which 
represents changes in the internal state. Thus, 


I 1 PSD 
w=—*t : d c A ~w+ tv 
Po 


D 

Only W enters the entropy production 


D 

TS = q+ W 


The second law of thermodynamics requires 
D 

W^O 


(34) 


(35) 
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We use as thermodvnamic state variables the elastic strain, 
£ 

represented by /£, the absolute temperature 7", and a number 
of other internal state variables ) that 

£ 

may be scalars and tensors of even order. The choice of f* 
and T as state variables is based on the fact that in pure ther- 
moelastic deformations, both quantities form a suitable set 
of thermodynamic state variables. The plastic strain and the 
total strain are unsuitable as state variables because, in 
general, they do not uniquely define the state of the 
material. A conflicting point of view has been expressed in 
Refs. 42-44. The remaining state variables are added for the 
sake of the description of the changes of the internal struc- 
ture of the material. 

The specific free energy (Helmholtz function) 0 given by 

<t> = u-Ts (37) 

must be a unique function of the thermodynamic state 

variables 

= (38) 

Since the elastic part of the deformation, according to our 

assumptions, does not depend on the plastic deformation, we 
may divide the free energy into two different components, as 

£ £ S 

0 = 0(/o T) +<t>( T, Av..,a£,.., # /4£f),...) (39) 

£ 

where the first component 0 refers to the elastic deformation 
S 

and the second 0 to the changes of the internal state. 

From Eqs. (31), (33-35), and (37) we derive 


deformation. Thus, we assume, in general 

D I 

W=C£ B D T c A d%>0 (43) 

where 

£ 

= (44) 

Equations (42) and (43) are the governing equations for 
nonisothermal, elastic-inelastic elementary processes. The 
specific free energy 0, which determines the nondissipated 
work of the thermomcchanical process and the quantity Qfg, 
which governs the entropy production, must be specified ac- 
cording to the material behavior. 

Eiastic-Viscoplastic Mode! 

Thermomechanical processes in elastic-viscoplastic bodies 
cannot be considered as a sequence of equilibrium states, 
even in the case of the elementary processes considered here. 
Elastic-viscoplastic deformations are associated with non- 
equilibrium states. One consequence of this fact is that we 
may get a continuation of a process without any change in 
the independent process variables. This occurs, for example, 
in the case of creep with constant stress and temperature or 
in the case of an adiabatic stress relaxation under constant 
strain. In such cases, the body moves from a nonequilibrium 
state to an equilibrium state. 

In order to establish the constitutive relations for elastic- 
viscoplastic bodies, which in the limiting case becoming 
elastic-inviscidly plastic, we adopt the usual assumption that 
the stresses, which produce the inelastic deformation, may be 
expressed as the sum of the so-called athermal or inviscid 
stresses, and the viscous overstresses r£ 

= (r£“ f£) (45) 


£ S 

0 = ~sr+ IV+ W 


(40) 


Also, we obtain from Eq. (39) 

£ £ E S 
'■0 Y 3(0 + 0) • 

* 3/2 Jc dT 
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By comparison of Eqs. (40) and (41), we may conclude that 


5= - 
S 


£ S 
3(0 + 0) 


dT 
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(42) 


For irreversible processes, this scheme of description has 
to be completed by some statements about the dependence of 
entropy production on the thermomechanical process. Under 
our assumption, we need deal only with entropy production 
by dissipated mechanical work, in connection with inelastic 


This assumption by no means detracts from the “unified” 
concept. The rate-independent limit of viscoplastic con- 
stitutive relation was recently discussed by Travnicek and 
Kratochvil. 45 Hence, the total work rate can be partitioned 
in the following way: 


. £ P V 1 £ i P i * P 

W= W+ W+ IV= — Ted a + — *cd C A + — T$d c A 
Po Po Po 

I 

W 


(46) 


The viscous part of the work is completely dissipated. Thus, 
we may write 

v Q, 

W=W (47) 

Regarding the plastic work, we have already stated that one 

part is used for changing the internal state and only the re- 
maining part can be considered to be dissipated. Therefore, 
we must write 

PSD p 

W= W+ W (48) 

So, we finally obtain 

£ S D D k 

iv= w+ ( 49 > 

D 

W 
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We have assumed that the changes of the internal state of 
the material can be regarded as a sequence of equilibrium 
states. Then, the specific energy is well defined in each state 
of the process and we may take the usual overall statement 
concerning the specific free energy. In so doing, however, we 

S 

must be aware of the fact that into the part W of the plastic 
P 

work rate W only the athermal stress enter, since only 
these stresses are involved in the plastic mechanism. For the 
same reason, we can introduce only the athermal stresses f£ 

D p 

into the statement concerning the dissipated plastic work W . 
On the other hand, we have to add the dissipated viscous 
D v D, 

work W to w in order to obtain the total rate of dissipation. 
The different mechanisms for determining the total dissipation 
and their coupling have been discussed by Perzyna. 46 

We now consider an example in which the specific free 
energy has the following form: 

EE S EE 

* = <*>( f A c. T) + <t)(T, k, a£)=4>(f£, T) + k +/( T) + ka£a% 

(50) 


as internal variables into the corresponding constitutive 
equations of the process description. 

The constitutive equations themselves are not yet deter- 
mined completely by Eqs. (50), (51), (55), and (56). These 
given only the restrictive frame for the formulation of these 
equations. We may derive a complete set of constitutive 
equations, which is compatible with this frame, by the fur- 
ther assumptions: 

1) The introduction of a yield condition of the form, 

P 

F~{t£- cp 0 hctt ) ( - cpo hct£ )-g 1 (W,T) =: 0 (57) 


where t£ denotes the deviator of the Kirchhoff stresses f£. 

2) The plastic deformation obeys the so-called normality 
rule, 


P . dF 


(58) 


3) The relations between the viscous stresses and the in- 
elastic deformation rate are of the form. 


P 

d A c 





(59) 


In this equation, h denotes a constant with the dimension of 
a specific energy like the variable k and the function /(T*). 
Furthermore, we assume that the dissipation is given by 

d p i p 

S(f£-cp 0 hc4)d$ 

Po 

D v 1 P 

W~ — (T^-r^)d^ (51) 

Po 


4) The quantities £ and c are constant. 

We can eliminate the athermal stresses f£ (which are not 
state variables) from the equations of evolution by consider- 
ing that the inelastic deformation can be expressed in two 
different ways. In one, the plastic mechanism is considered 
and the viscous mechanism in the second. From Eq. (57), we 
then obtain 

P 

d£ = 2 \(f£ — c/> 0 hctc ) (60) 


where £ < 1 and c denotes constant numbers. This leads to 

D EL D v P PI 

w + (£- l)If r -£c/ta£c£ + W (52) 


Hence, we obtain 

SID P 

W=W-W=(\~-t)iv+tchcx£dZ (53) 

On the other hand, from Eqs. (42) and (50) we have 

s v 

W=k + 2hcc£a% (54) 

Equations (53) and (54) are compatible, for instance, if we 
put 


P 


k={\-l)W 

(55) 

V P 

<*A = ‘/iCStfS 

(56) 


P 

From Eq. (55) it follows that, in our case, the plastic work If 
is equivalent to the thermodynamic state variable k . This is 
still true if we take £ as a function of k . But it does not hold 
in the general case when £ also depends on the other state 
variables T and a £. Equation (56) shows that only in a very 
special case, a very unrealistic one, the state variable is 
equivalent to the plastic deformation. 

From the thermodynamical considerations, it follows then 
that we may introduce the quantities k and a£, defined by 

P 

Eqs. (55) and (56) or any other equivalent set (W,cp 0 hot£)\ 


while from Eq. (59), we have 

^,J_(^-^) 

2ri 

= (/£ — ( fc ) ) (61) 

p 

By comparing these equations for we get 

\ JL £ ^ (t£~-cp 0 ha£) (t^—cpphot^ * - 1 j (62) 

Following the course of the process in each state, the internal 
P P 

parameters W and ct£ and, therefore, also k ? = k 2 (W, «£) are 
known. Thus, we may calculate X from Eq. (62) and then all 

P 

the other needed quantities such as and 

Discussion 

Many thermodynamic considerations of nonisothermal, 
elastic-viscoplastic deformations refer essentially to the 
general fundamentals that must be observed in describing 
such phenomena as thermomechanical processes and then 
discuss what particular restrictions follow from the second 
law of thermodynamics. Only a few papers attempt to 
describe completely such processes by state equations. Most 
of these papers introduce plastic strains as thermodynamic 
state variables. But one may conclude from the consideration 
of the phenomena in the crystal lattice (dislocations, for ex- 
ample, that have completely passed through the crystal pro- 
duce plastic strains but no changes of state) as well as from 
phenomenological observations (different states of hardening 
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can belong to the same plastic strains) that plastic strains in 
general cannot be regarded as state variables. Furthermore, 
all these papers consider the plastic work as completely 
dissipated. However, this is in contradiction with experimen- 
tal results, from which it emerges that one part of plastic 
work is used for producing states of residual stresses in the 
lattice, which, when phenomenologically considered, cause 
hardening. 

The results in work presented here can be extended to 
more complex constitutive equations by introducing more in- 
ternal parameters or state variables. We may extend our ap- 
proach to more general, anisotropic hardening materials by 
assuming [see Eq. (50)], for example, that 

EE S 

<i> = C/o T) +o( T,k,a^ f AcD) 


E E 

= <t> (. T) + k +/( T) +A£ B D a c A a D B (63) 

Also, it may be more advantageous to replace the assump- 
tion in Eq. (58) for the plastic deformation rate by 

P dF v 

+ (64) 

Off 

This form of this model appears to be more suitable for 
representing some experimental results in w'hich second-order 
effects and some deviations from the normality rule have 
been observed. Sometimes, the normality rule is considered 
as a fundamental law based on an entropy production princi- 
ple. But we should keep in mind that, since not all of the 
plastic work is dissipated, we cannot expect the total plastic 
deformation rate to obey the theory of plastic potential even 
though the mentioned principles of entropy production are 
correct. 
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Nonisothermal Elastoviscoplastic Snap-Through and Creep 

Buckling of Shallow Arches 


R. Riff* 

Georgia Institute of Technology , Atlanta, Georgia 

\ 


The problem of buckling of shallow arches under transient thermomechan.cal load! « • '”**'* ’ ™ 

analvsis is based on nonlinear geometric and constitutive relations and is expressed in a rate fonn. Tbe maten 
constitutive equations are capable of reproducing all nonisothermal. elastoviscoplastic characteristics. The 
Son scheme is capable of predicting response that includes pre- and postbuckling with creep and p a tic 
effects. The solution procedure is demonstrated through several examples that include both creep and snap- 
through behavior. 


Introduction 

E LASTIC snap-through of low arches under quasistatic 
loads has been the subject of several investigations. The 
significance of the problem, insofar as it illustrates- certain 
important features in more complicated buckling problems of 
piastes and shells, was pointed out by Marguerre, 1 who con- 
structed a simplined mechanical model to demonstrate these 
features. Timoshenko* obtained an approximate solution to 
the problem of a low arch under a uniformly distributed load. - 
Biezeno 5 considered the problem of a low-parabolic arch 
loaded laterally at the midpoint with a concentrated load. He 
also investicated snap-through buckling of a shallow spherical 
cap, pinned* along its circular boundary, under the action of a 
concentrated load applied along the axis of rotational symme- 
try* He presented his approximate solutions by means of load- 
deflection curves and equations from which the critical load 
could be calculated. 

In 1952, Fung and Kaplan 4 investigated the problem of low- 
pinned arches of various initial shapes and spatial distributions 
of the lateral load. Their results show that a very shallow arch 
snap through symmetrically, whereas a higher arch buckles 
asvmmetrically. They also ran a limited number of experimen- 
tal tests, and 'their experimental data are in good agreement 
with their theoretical results. About the same time, Hoff and 
Bruce 5 , in investigating the possibility of snap-through buck- 
line of a low-pinned arch with a half-sine-wave initial shape 
under a half-sine-wave distributed dynamic load (suddenly ap- 
plied with infinite duration), obtained results for the qua- 
sistatic load case that are identical to those obtained by Fung 
and Kaplan for the same problem. 

In 1962, Gjelsvik and Bodner* obtained an approximate 
solution to the problem of a low-circular arch with a concen- 
trated load at the midpoint of the arch and clamped boundary 
conditions. They also reported on experimental results. 
Schreyer and Masur obtained an exact solution to their prob- 
lem (and for the load case uniform pressure), and they showed 
that for the concentrated load case, the arch snaps through 
symmetrically regardless of the value of the rise parameter. 
Masur and Lo* presented a general discussion of the behavior 
of the shallow circular arch regarding buckling* postbuckling, 
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and imperfection sensitivity. Snapping of low-pinned arches 
resting on an elastic foundation has been investigated by Sim- 
iises. 9 " This system exhibits all forms of experimentally ob- 
served buckling phenomena (smooth and violent) and of theo- 
retically predicted responses (limit point, bifurcation with 
stable branching, and bifurcation with unstable branching), 
and it is presented with sufficient detail in Ref. 10. Experimen- 
tal results have also been reported by Roorda. 11 

The effects of inelastic material behavior found theit way 
into the literature since the 1960’s. Onat and Shu :: considered 
the behavior to be one of rigid-perfectly plastic. Fromciosi, 
Augusti, and Sparacio 13 discussed the collapse of arches under 
repeated loads with inelastic material behavior. Studies of in- 
elastic snap-through buckling of shallow arches also^were re- 
ported by Lee and Murphy. 14 In addition, August! 15 ' investi- 
gated plastic buckling of a model of a .three-hinged aren in 
1968, and a more complete analysis of the same model was 
provided by Batterman 16 in 1971. Finally, the reader who is 
interested in the ultimate strength of parabolic steel arches 
with b ra cine system is referred to Komatsu,* who considers 
inelastic in-plane and out-of-plane instabilities and provides 
design formula for each case. 

Creep buckling of shallow arches has been investigated by 
Huang and Nachbar, 1 * Krajcinovic, 19 and Botros and Bi- 
enek. :o The elastic response of arches under sudden (dynamic) 
application of the external loads has been reported by Hoff 
and Bruce. 5 Hsu. :i ~ and Lock. B For a more complete bibliog- 
raphy see Ref. 24. As far as the authors knew, no work has 
been reported on the nonisothermal elastoviscoplastic behav- 
ior of shallow arches. The purpose of this paper is to demon- 
strate ihe effect of highly nonlinear material behavior on the 
snap through and creep buckling of shallow arches. 

Elastothermoviscoplastic Constitutive Relations 
The prediction of buckling loads and postbuckling behavior of 
structural components, like shallow arches, made of a realistic 
metallic material and subjected to nonisothermal thermome- 
chanical loads has increased in importance in recent years. 

Under this kind of severe loading conditions, the structural 
behavior is highly nonlinear due to the combined action of 
geometrical and material nonlinearities. On one side, finite 
deformation in a stressed structure introduces nonlinear geo- 
metric effects. On the other side, physical nonlir.eariiies arise 
even in small strain regimes, whereby inelastic phenomena 
plav a particularly important role. From a theoretical stand- 
point, nonlinear constitutive equations should be applied only 
in connection with nonlinear transformation mta'U'r: ?.".r'* • 
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ing both deformation and rotations). However, in almost all of 
f the works in this area, :5 the two identified sources of nonlin- 
•earictes are always separated. The separation yields* at one end 
of the spectrum, problems of large response, whereas at the 
other end problems of viscous and/ or nonisothermal behavior 
in the presence of small strain. 

The classical theories in which :he material response is char- 
acterized as a combination of distinct elastic* thermal, time-in- 
dependent inelastic (plastic) and time -dependent, inelastic 
(creep), deformation components cannot explain some phe- 
nomena that can be observed in complex the rmo mechanical 
loading histories. This is particularly true when high- tempera- 
ture nonisothermal processes must be taken into account. 
There is a sizeable body of literature 2526 an phenomenological 
constitutive equations for the rate and temperature-dependent 
piastic deformation behavior of metallic materials. However* 
almost all of these new “unified” theories are based on small 
strain theories* and several suffer from some thermodynamic 
inconsistencies. 

In a previous paper* 2 * the authors have presented a complete 
sec of constitutive relations for norisothermaf* targe strain* 
eiascoviscopiasExe behavior of metals. It was shown, there 27 that 
the metric censor in the conveeted (material) coordinate system 
can be linearly decomposed into elastic and (visca) plastic 
parts. So a yield, function, was assumed, which is: dependent on 
the race of change of stress on the metric* on the temperature, 
and on a set of internal variables. Moreover, a hypoeiastic 
law was chosen to describe the rhermoeiasctc part of the 
deformation.. 

A time- and temperature-dependent vtscopl'asticity model 
‘ was formulated in this conveeted material system to account 
for finite strains and rotations. The history and temperature 
dependence were incorporated through the introduction of 
internal variables. The choice of these variables, and their evo- 
lution were motivated by thermodynamic considerations. 

Jhe nonisothermal elastov isc o pi as t ic deformation process 
was described completely by “thermodynamic state” equa- 
tions. Most investigators 1 - 26 (in the area of viscoplasticity) 
employ piastic strains: as state, variables. The authors* previous 
study- 7 shows, chat, in general, use of plastic strains as state 
variables may lead to inconsistencies with regard to- thermody- 
namic considerations. Furthermore* the approach and formu- 
lation employed in previous works lead m the condition that 
ail of the plastic work is completely dissipated. This*, however* 
is in contradiction with experimental evidence* from which it 
emerges that part of the plastic work is used far producing 
residual stresses in. the. take that, when phenomenologically 
considered, causes hardening. Both, limiarions were excluded 
from, this 27 formulation. Accuracy of the formulation was 
checked on a wide raqge. of examples. :s 

The constitutive relation will be rephrased here in some 
different form. For brevity we compile only some notations 
and fundamental relations that are used Ins the formulation of 
the incensed constitutive law. For details see Refs. IT and 2SL 

Concerning; the formulation: of constitudve laws* k is advan- 
tageous: to use a material (conveeted) coordinate system. The 
transformation from the under formed state (metric £*) to the 
deformed state g 4k can be represented by the tensor 


and u iik are the material velocities gradients. Hence, 

d :k m Vi(v Lk +■ u k ',) (4) 

The expression 

h = (/)'* - 4/* - difi * sym(CO^J . (5) 

represents; the symmetric part of (/)'* or the covariant deriva- 
tive with respect to time, also called the conveeted derivation, 
which is due to Zaremca and Jaumann. 

The total deformation can be decomposed according to 

* . torn 

fl=i im 3«ris*=xn (6) 

where the superscript (*) relates to a fictitious configuration 
defined by a fictitious reversible process with frozen internal 
variables. The decomposition of Eq. (6) leads to an additive 
decomposition of the deformation rate 

(n (t) 

4 - <4 -r d i (7) 

W (0. 

d\ is related to the reversible deformations, whereas d k de- 
notes the remaining part of the deformation rate. 

For the description of the stress state, we introduce the 
Kkchoff stress tensor si, which is connected with the real 
Cauchy stress tensor a : k by the relation 

si = (p/afai ( 8 ) 

Assuming that the elastic behavior is not affected essentially 
by inelastic deformations, the following hypoeiastic incremen- 
tal law was chosen 17 : 



where c k is the weighted stress deviate r, C the shear modulus, 
K the bulk modulus, and a the coefficient of thermal expan- 
sion. 

The following constitutive relations were established 27 for 
the inelastic behavior. Yield condition: 

F ~itl- chtiXj? - T) - * 2 >0 : (10) 

Accompanying equilibrium state: 

f - cfatixi? - chsti - n = /- - * 2 = o (i i) 

Evolution law for inelastic deformations:. 

m 

d\ = 2X(7; - (12) 

with 


Jit — 3 r ZrK 

or 

(la) 

tjm-cpgmtf -csgisn „ 

X = 4, (N l) 

(13) 


(lb) 

7* = j + + cpgPic 

(14) 

The total deformation rate is defined by 


Evolution laws for the internal variables: 


4 * '/ir'g fk = = Hi r'TUfr* - 

(2) 

A = t t%d\ 


Here ( * ) denotes time material derivative* The rate of change 

(15) 

of the metric tensor is given; by 


V <i)> 


&<k, - V,.k - f*., 


4 - 

(16) 
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If 


F=0and IJ^ + l? f>0 

(17) 

«•> 


di = d i 

(18) 

d\ = 0 and d [ = 2X(r; - cpgfr) 

(19) 

x-g -p ZUi-cpgfalf — T 

(20) 


If 

BF v . BF ■ 

F = Q and — s‘ K ^ — 7" < 0 (21) 

05a 0 1 

or 


F< 0 

(22) 


II 

(23a) 

A = 0 

(23b) 

! 

V. ■ 

=0 

(23c) 


Within the developed frame, the elastoviscoplastic behavior 
is governed bv the scalar material functions c(s l k ,T,A, ft), 
Ar 2 (k,D, *<5'>,ft), £M,A,ft), and ^,7, ft). These ma- 
terial functions can be determined from a series of monotonic 
and cyclic processes with proportional and nonproportional 
paths at different temperature levels. 25 

General Formulation and Solution Schemes 
The rate form of the constitutive equations suggests that a 
rate approach be taken toward the entire problem so that flow 
is viewed as history-dependent process rather than an event. 
For this purpose, a complete true ab-initio rate theory of kine- 
matics and kinetics for continuum and curved thin structures, 
without any restriction on the magnitude of the transforma- 
tion, was presented in Ref. 28. It is implemented with respect 
to a body-fixed system of convected coordinates, and it is valid 
for finite strains and finite rotations. The time dependence and 
large strain behavior are incorporated through the introduc- 
tion of the time rates of change of the metric and of the 

curvature. 

The constitutive law may be applied to the conservation of 
momentum via an appropriate variational principle. The prin- 
ciple of virtual power (or of virtual velocities) reads 

dV- \^vpf-bvj dV- \ a »T'&vj dA = 0 >(24) 


where 6v y are the virtual velocities, p the body forces per unit 
mass, and vT J the surface tractions. Total differentiation of 
Eq. (24) yields 



At any instant Eq. (25) must be satisfied. The virtual veloc- 
ity and its time derivative are then independent. Moreover, the 
last three terms of Eq. (25) are equivalent to Eq. (24). Hence, 
the principle of the rate of virtual power may be obtained in its 
concise form. For further classifiations, the total derivative of 
the stress components will be represented by the Jauman 
derivative, and the following integrals are defined by 

I t = l <j" ovjj d V (26) 

J*' 

h = [ - a^Yovj, dV (27) 

J v 

ft 

/, = \yi > ^bvj, i dV ( 28 ) 

Then, substitution in Eqs. (23) yields the final form of the- 
principle of the rate of virtual power: 

f dp f dP 

/ = /,r/ c t/ f * | hvj dV + V— ovj dA (29) 

The auasiiinear nature of the principle of the rate of virtual 
power suggests the adoption of an incremental approach to 
numerical integration with respect to time. The availability of 
the field formulation provides assurance of the completeness 
of the incremental equations and allows the use of any conve- 
nient procedure for spatial integration over the domain V. In 
the present instance, the choice has been made in favor of a 
simple first-order expansion in time for the construction of 
incremental solutions from the results of finite-element spatial 
integration of the governing equations. 

The procedure employed permits the rates of the field for- 
mulation to be interpreted as increments in the numerical solu- 
tion. This is particularly convenient for the construction of 
incremental boundary condition histories. 

The finite-element method for spatial discretization has 
been well documented (see, e.g., the books by Zienkiewicz 29 or' 
Oden 50 ) ana will not be detailed here. The algebraic counter- 
pan of Eg. (29) After the finite-element discretization (for 
detail see Ref. 28) is the well-known stiffness expression 

[A'ltn «{/>)- IF} (30) 

with the tangent stiffness matrix [A'], the vector of the incre- 
mental velocities (I'), and the vector out-of-baiance force 
rates, external force rates [P] minus internal force rates [F j. 
The homogenous case of Eq. (30) indicates either the non- 
uniqueness of the equilibrium path at a stable or unstable 
bifurcation point or the unique but unstable situation at a limit 
point. Hence, this criterion may be evaluated by a determinant 
check or supplementary eigenvalue analysis for the load 
parameter parallel to the loading process. 

Even under the condition of static externa! loads and slowly 
growing creep effects, the presence of snap-through buckling 
makes the inertia effects significant. In dynamic analyses, the 
applied body forces include inertia forces. Assuming that the 
mass of the body considered is preserved, the mass matrix can 
be evaluated prior to the time integration using the initial 
configuration. 

Finite-element solution of any boundary-value problem in- 
volves the solution of the equilibrium equation (global) to- 
gether with the constitutive equation (local). Both equations 
are solved simultaneously in a step by step manner. The incre- 
mental form of the global and local equations can be achieved 
by taking the integration over the incremental time step 
A / = l Jm , - t,. The rectangular rule has been applied to execute 
the resulting time integration. 

Cleariy, the numerical solution involves iteration. A simpli- 
fied version : of Riks Wemr r.er cor.::cr.t-a r ch-!er:g:h method 
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Fig. 1 damped circular arch. 



Fig. 1 Paraiinear isoparametric element. 

has been used. This iteration procedure, which is a generaliza- 
tion of the displacement control method, also allows to trace 
the nonlinear response beyond bifurcation points. In contrast 
to the conventional Newton- Raphson techniques, the iteration 
of the method takes place in the velocity and load rate space. 
The load step of the first solution in each increment is limited 
by controlling the length ds of the tangent. Either the length is 
kept constant in each step, or it is adapted to the characteristics 
of the solution. In each step the triangular-sized stiffness ma- 
trix has to be checked tor negative diagonal terms, indicating 
that a critical point is reached. 

Shallow Circular Clamped Arch 

The theory and computational procedures described in the 
proceeding sections have been applied to the creep buckling 
analysis of a shallow circular clamped arch. The problem of 
the clamped arch, besides being of some practical interest, 
contains a number of similarities to that of the uniformly 
loaded spherical cap. The trend of results of the arch problem 
serves as a good indicator to the behavior of the latter. 

The shallow circular clamped arch subjected to a single cen- 
tral concentrated load, as shown in Fig. 1, is analyzed. The 
material chosen for the numerical experimentation is the car- 
bon steel C-45 (DIN 1720) with £ = 10 T psi, v - 0.3, and 
a , ~ 2.7 x 10 4 psi at room temperature. The viscoplastic prop- 
erties (the scalar material functions) were obtained in Ref. 28. 

The analysis is performed with the aid of 24 paraiinear 
isoparametric elements (Fig. 2). The paraiinear isoparametric 
element is intended for the analysis of oriented structures 
where the geometry is such that the thickness is small com- 
pared to other dimensions. The characteristics of the element 
are defined by the geometry and interpolation functions, 
which are linear in the thickness direction and parabolic in the 
longitudinal direction (see Fig. 2). Consequently, the element 
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Fig, 3 Arch response. 
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Ftg. 4 Deflection time history. 

allows for shear strain energy, since normals to a midsurface 
before deformation remain straight but not necessarily nor- 
mal to the midsurface after deformation. 

The elastic behavior, corresponding to both axisymmecric 
and asymmetric response, is shown on Fig. 3. These curves are 
in complete agreement with those produced by Gjeisvik and 
Bodner, 6 only because the Young's modulus and Poisson’s 
ratio values used are virtually the same (carbon steel C-45 here 
and 2024-T4 aluminum alloy in Ref. 6). Note that these elastic 
response curves are hypothetical for our material but true for 
the 2024-T4 alloy. The true behavior for our material is elasto- 
viscoplastic, and it is labeled as such on Fig. 3. Note that this 
curve represents quasistatic (steady-state) eiastoviscoplastic re- 
sponse. as described by the chosen constitutive law. According 
to this, snapping occurs at a load of 26.20 lb, primarily be- 
cause of the low-vield strength. Then, the postlimit point be- 
havior seems to be primarily driven by viscoplastic behavior. 

It is expected here that if loads up to approximately 14 lb are 
reached quasistatically and left applied for a long time, the 
primary response will be creep, and the critical creep collapse 
time will be extremely large. On the other hand, for loads 
between 14 and 26.2 lb (especially for the higher range), the 
behavior will be a combination of creep and snap-through 
buckling. This is best demonstrated by the curve on Fig. 4. The 
applied load is reached quasistatically in 13 min, and then it is 
kept constant. The curve of Fig. 4 depicts the behavior of the 
arch in terms of midpoint deflection vs time. Creep, initially, 
is very slow; then snap-through takes place in 32 min, curve 
BC, and then the creep behavior continues until a critical time 
to creep (creep buckling occurs) is reached after a total time 97 
min. Note that for this loading condition, the critical time to 
creep in 97 min. Creep buckling and critical time to creep are 
based on the phenomenon that the deflection increased very 
rapidly. For loads higher than 26.2 lb, it is expected that snap- 
ping will occur early during quasistatic loading, and then the 
creep behavior will be similar to that shown on Fig. 4, past 
point c. ' . 







Fig. 6 Influence of temperature raise. 


The next example considers the influence of cyclic loading 
on the response. Figure 5 illustrates the load deflection behav- 
ior of the arch under a cyclically applied external load. The 
load is increased, quasistaticallv, from zero to 26 lb in 5 min, 
then it is held constant for 2.5 min. after that it is reduced to 
20 lbs and held constant for 50 min, then raised to 25.5 lb for 
2.5 min, and finally it is reduced back to 20 lb held constant. 

The steady-state response under this type of loading exhibits 
several relative maxima points, which may imply that snapping 
is imminent shortly after the load reaches the value of 26 lb 
(between points A and B on Fig. 5). The dashed curve corre- 
sponds to the hypothetical elastic static response, and it is only 
shown for comparison purposes. 

The last example presented in Fig. 6 considers the influence 
of temperature on the arch behavior. The loading history is the 
same on the one shown on Fig. 4. The curve corresponding to 
r- 50*F was discussed previously i Fig. 4), and it is used here 
as a basis for comparison. When the temperature is raised to 
200*F (after this, the loading is applied), the time to snap is 
reduced to 26 min, whereas the critical creep collapse time is 
not appreciably affected. On the other hand, at the highest 
temperature T - 1000*F for which results are shown. The crit- 
ical creep collapse time is reduced to 62 min, and the steady- 
state response does not show a clear snap-through behavior. 
Different values of a were used for the different temperature 
in the elastic range. 


Discussion 

As noted earlier, the main thrust of this work has been to 
demonstrate the effect of highly nonlinear material behavior 
on the snap-through and creep buckling of shallow arches. It 
is evident that in the presence of both elastic and viscoplastic 
deformation the process of buckling assumes an entirely new 
character. Buckling develops as a time-temperature-dependent 


deformation process under constant or variable loads of mag- 
nitudes smaller than the elastic critical values. In arenes under 
loads below the critical values the structure initially deforms 
quasistatically, with the thermoviscous terms manifesting 
themselves in the form of increasing displacement under, say, 
a constant load. When the magnitude of the displacements 
reaches a certain threshold state, the arch snaps dynamically 
into the postbuckling configuration and then continues qua- 
sistatic deformation again. 

The material constitutive relation has been proven to be 
capable of reproducing the main characteristics ot viscoplastic 
deformations. The modified Riks/Wempner iteration scheme 
has been found to be a versatile technique in the pre- and 
postcritical range. 

The influence of thermomechanical coupling can become 
very large in stability problems. Such processes are always 
connected with a rapid grow-th of inhomogeneity of the state. 
Thorough investigation of such problems, however, must be 
performed with the necessary detail. 
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NON-ISOTHERM AL BUCKLING BEHAVIOR OF VISCOPLASTIC SHELL 

STRUCTURES 


Richard Riff 

School of Aerospace Engineering 
Georgia Institute of Technology _ 

Atlanta, Georgia 30332 

Abstract 

The buckling behavior of thin metallic shell structures under transient thermomechani- 
cal loads is investigated. The analysis is based on nonlinear geometric and constitutive 
relations, and all the field equations are expressed in a rate form. The employed consti- 
tutive equations are thermodynamically consistent and they are capable of capturing all 
non-iso thermal, elasto-viscoplastic characteristics of the response. The solution scheme 
is capable of predicting response which includes pre- and post— buckling with creep and 
plastic effects. The solution procedure is demonstrated through several examples which 
include both creep and snap-through behavior. 

Introduction 

The prediction of inelastic behavior of metallic materials at elevated temperatures has increased 
in importance in recent years. The operating conditions within the hot section of a rocket motor 
or a modern gas turbine engine present an extremely harsh thermomechanical environment. Large 
thermal transients are induced each time the engine is started or shut down. Additional thermal 
transient from an elevated ambient, occur whenever the engine power level is adjusted to meet 
flight requirements. The structural elements employed to construct such hot sections, as well as 
any engine components located therein, must be capable of withstanding such extreme conditions. 
Failure of a component would, due to the critical nature of the hot section, lead to and immediate 
and catastrophic loss in power and thus cannot be tolerated. Consequently, assuring satisfactory 
long term perfor- mance for such components is a major concern for the designer. 

The problem of inelastic analysis of shell structures has been investigated recently by Kojic’ 
and Bathe 1 . They used the "effective-stress-function” algorithm to compute plastic and creep 
effects on the behavior of shell like structures. The effects of inelastic material behavior on stability 
of shells found their way into the literature since the late 1970’s. The paper by Miyazaki and 
Ando 2 deals with creep buckling of perfect spherical shells subjected to pressure loading and 
considers only the effects of steady-state creep. Xirochakis and Jones 3 have studied axisymmetric 
and bifurcation creep buckling of externally pressurized spherical shells under the condition of 
secondary creep only. Botros and Bienek 4 presented a numerical treatment of the creep buckling 
of these configurations. Their work includes the effects of elastic strain, primary and secondary 
creep strains and creep recovery. The influence of temperature and viscous effects on dynamic 
buckling of shells has been considered by Wojewodzki and Bukowski 5 ’ 6 . The book by Owen and 
Hinton 7 gives a list of references for the applications of finite element methods to the problem of 
creep buckling of structures. 

Assistant Professor, Member AIAA 
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As far as the author know no work has been reported on the non-isothermal buckling behavior of 
elasto-viscoplastic shell structures.The purpose of this paper is to demonstrate the effect of highly 
nonlinear material behavior on the snap through and creep buckling of shells. 


Elasto-Thermo— Viscoplastic Constitutive Relations 


In a previous works 8 -“, following the ideas of Lechmann 12 - 13 the authors have presented a 
complete set of constitutive relations for nonisothermal, large strain, elasto-viscoplastic behavior of 
metals. It was shown there 8 that the metric tensor in the convected (material) coordinate system 
can be linearly decomposed into elastic and (visco) plastic parts. So a yield function was assumed, 
which is dependent on the rate of change of stress on the metric, on the temperature and on a set 
of internal variables. Moreover, a hypoelastic law was chosen to describe the thermo-elastic part 

of the deformation. 


A time and temperature dependent viscoplasticity model was formulated in this convected 
material system to account for finite strains and rotations. The history and temperature dependence 
were incorporated through the introduction of internal variables. The choice of these variables, as 
well as their evolution, was motivated by thermodynamic considerations. 

The constitutive relation will be rephrased here in some different form. For brevity we compile 
only some notations and fundamental relations which are used in the formulation of the intended 
constitutive law. For details, see Refs. 8 and 11. 

Concerning the formulation of constitutive laws it is advantageous to use a material (convected) 
coordinate system. The transformation from the underformed state (metric <7,*) to the deformed 
state ( gik ) can be represented by the tensor: 


fi=9 ir 9rk or (r 1 )i = 3‘ > 5r* W 

The total deformation rate is defined by 

4 = = |( r'YAm = -|(r l )U < 2 ) 

here f) denotes time material derivative. The rate of change of the metric tensor is given by 

9ik = v i,k + v k,i 

and v^jfc are the material velocities gradients. Hence, 

<*«* = |(vf,fc + Vj|.,,) ( 4 ) 

The expression 

represents the symmetric part of (/)** or the covariant derivative with respect to time, also called 
the convected derivation, which is due to Zaremba and Jaumann. 
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The total deformation can be decomposed according to 


? * 

fl =9 ,m 9mr 9, k =f'rf r k 


( 6 ) 


where the superscript ( ) relates to a fictitious configuration defined by a fictitious reversible process 
with frozen internal variables. The decomposition of Eq. (6) leads to an additive decomposition of 
the deformation rate 

If I HJ 

(7) 


(r) (0 

4 =4 + 4 


cfjj. is related to the reversible deformations, whereas denotes the remaining part of the deforma- 
tion rate. 

For the description of the stress state, we introduce the Kirchoff stress tensor sj., which is 
connected with the real Cauchy stress tensor a\, by the relation: 


A _ Pi 

*k - ~*k 


( 8 ) 


Assuming that the elastic behavior is- not affected essentially by inelastic deformations, the 
following hypoelastic incremental law was chosen 8 


to I v 1 . . 

4= i +■ 


with 


G 

K 

a 


weighted stress deviator 
shear modulus 
bulk modulus 

coefficient of thermal expansion 


The following constitutive relations were established 8 for the inelastic behavior, 
yield condition: 

F = (4 - cP - cP gfif) - k 2 {A,T) = f 2 - k 2 > 0 

accompanying equilibrium state: 

F = (4 - c P - C P gfi) - k 2 (A,T) = f - k 2 = 0 

evolution law for inelastic deformations: 

(0 


4= 2A(? fc ~ c P 90k) 


with 


and 


x = 1 l Z l \ - g j 


tl = 


1 + 4r?A 


(ti -cP gp{) + cPgp{ 


(9) 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 
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evolution laws for the internal variables? 


if 

then 

with 

if 

or 

then 


^ — a 4 4 

P 

(15) 

V (*) 

4=e4 

(16) 

&F v - BF ■ 

f = 0 and 5^*i+ay r>0 

(17) 

II 

(18) 

£«} V . ... o . 

4= 0 and 4= 2A(4 — c P } 

(19) 


(20) 

dF * BF • 

f|4 + ^T<0 

(21) 

F < Q 

(22) 

O 

II II 

(23) 


v 

fii = 0 


Within the developed frame the eiasto-viscoplastic behavior is governed by the scalar mate- 
rial functions c(4,r,A4>* •** n(AT,&). These material 

functions can be determined from a series of monotonie and cyclic processes with proportional and 
nonproportional paths at different temperature levels 

Formulation and Solution Schemes 

The rate form of the constitutive equations suggests that a rate approach be taken toward the 
entire problem so that flow is viewed as history dependent process rather than an event. For this 
purpose, a complete true ab-initio rate theory of kinematics and kinetics for continuum and curved 
thin structures, without any restriction on the magnitude of the transformation was presented in 
Ref. 11. It is implemented with respect to a body- fixed system of convected coordinates, and 
it is valid for finite strains and finite rotations. The time dependence and large strain behavior 
are incorporated through the introduction of the time rates of change of the metric and of the 
curvature. 



A nonlinear, thermodynamic theory of shells was derived in ref. 11 from three dimensional 
continuum mechanics in a natural and comprehensive way. All the approximations had been thrown 
into a postulated two dimensional form of the first law of thermodynamics. 

The derived shell theory, in the least restricted form, before any simplifying assumptions are 
imposed, has the following desirable features: 

(a) The two-dimensional, impulse-integral form of the equations of motions and the Second Law of 

Thermodynamics (Clausius- Duhem inequality) for a shell follow naturally and exactly from 
their three-dimensional counterparts. 

(b) Unique and concrete definitions of shell variables such els stress resultants and couples, rate of 

deformation, spin and entropy resultants can be obtained in terms of weighted integrals of 
the three-dimensional quantities through the thickness. 

(c) There are no series expansions in the thickness direction. 

(d) There is no need for making use of the Kirchhoff Hypotheses in the kinematics. 

(e) All approximations can be postponed until the descretization process of the integral forms of 

the First Law of Thermo- dynamics. 

(f ) A by-product of the descent from three-dimensional theory is that the two-dimensional tem- 

perature field (that emerges) is not a through-the-thickness average, but a true point by point 
distribution. This is contrary to what one finds in the literature concerning thermal stresses 
in the shell. 

The obtained complete rate field equations consist of the principles of the rate of the virtual 
power and the rate of conservation of energy, of the constitutive relations, and of boundary and 
initial conditions. These equations provide a sound basis for the formulation of the adopted finite 
element solution procedures. 

The quasi-linear nature of the principle of the rate of virtual power suggests the adoption of an 
incremental approach to numerical integration with respect to time. The availability of the field 
formulation provides assurance of the completeness of the incremental equations and allows the use 
of any convenient procedure for spatial integration over the domain V. In the present instance, the 
choice has been made in favor of a simple expansion in time for the construction of incremental 
solutions from the results of finite element spatial integration of the governing equations. 

The procedure employed permits the rates of the field formulation to be interpreted as incre- 
ments in the numerical solution. This is particularly convenient for the construction of incremental 
boundary condition histories. 

To develop geometrically nonlinear, doubly curved finite shell elements the basic equations of 
nonlinear shell theories have to be transferred into the finite element model. As these equations 
in general are written in tensor notation, their implementation into the finite element matrix 
formulation re- quires considerable effort. The nonlinear element matrices are directly derived 
from the incrementally formulated nonlinear shell equations, by using a tensor-oriented procedure. 
For this formulation, we use the ” natural” degrees of freedom per mid-surface shell node: three 
incremental velocities and the rates of rotations about the material coordinates in a mixed form. 
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Finite element solution of any boundary-value problem involves the solution of the equilibrium 
equations (global) together with the constitutive equations (local). Both sets of equations are solved 
simultaneously in a step by step manner. The incremental form of the global and local equations 
can be achieved by taking the integration over the incremental time step At = t J+ 1 — tj. The 
rectangular rule has been applied to execute the resulting time integration. 

Clearly, the numerical solution involves iteration. A simplified version of the Riks- Wempner 14 
constant-arc-length method has been utilized. This iteration procedure which is a genera- lization 
of the displacement control method also allows to trace the non-linear response beyond bifurcation 
points. In contrast to the conventional Newton-Raphson techniques, the iteration of the method 
takes place in the velocity and load rate space. The load step of the first solution in each increment 
is limited by controlling the length ds of the tangent. Either the length is kept constant in each 
step or it is adapted to the characteris- tics of the solution. In each step the triangular-size stiffness 
matrix has to be checked for negative diagonal terms, indicating that a critical point is reached. 

Applications 

Two different material representing different sensitivity to creep and high temperature Were 
chosen for the numerical experimitations. 

The first is the carbon steel C-45 (DIN 1720) with E = 10 7 psi, u = 0.3 and cr v = 2.7 x 10 4 psi 
at room temperature. The viscoplastic properties (the scalar, material functions) were obtained in 
Ref. 11. The C-45 has low resistance to temperature and high strain rate sensitivity. The second is 
Hastelloy-X a nickel-base alloy used for burner-liner parts, turbine-exhaust weldments, afterburner 
parts, and other parts requiring high resistance to temperature effects. Material properties used 
for Hastelloy-X were based on data from Refs. 15 and 16. 

The first example is of hinged spherical cap made of C-45 carbon steel and loaded by a con- 
centrated force at the appex (Fig. 1). The hypothetical elastic response is shown in Fig. 2 by the 
dash line. The true behavior for our material is elasto-viscoplastic, the full curve in Fig. 2. Note 
that this curve represents quasi-static (steady state) elasto-viscoplastic response, as described by 
the constitutive law. According to this, snapping occurs at a load of 32.3 lbs, primarily because 
of the low yield strenght. Then, the post-limit point behavior seems to be driven by viscoplastic 
behavior. 

It is expected here that for loads between 4 lbs and 32.3 lbs (especially for the higher range) 
the behavior will be a combination of creep and snap-through buckling. This is demonstrated by 
Fig. 3 for different temperatures. The applied load of 25 lbs is reached quasi-statically and then it 
is kept constant. Note the different creep buckling behavior and the differnt critical time to creep 
for the different temperatures. 

The second example considers a clamped spherical cap subjected to uniform pressure and made 
from Hastelloy-X (Fig. 4). Fig. 5 describes the different load deflection curves of the center point 
for different temperatures for this example. The critical creep collapse time is reduced by 45 

In the next example (Fig. 7) the central deflection time history and the influence of temperature 
change of a thin, imperfect, cylind- rical shell panel made of carbon steel C-45 is shown. The panel 
(Fig. 6) is simply supported on all sides, and subjected to inplane loads along the generator, The 
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Fig. 1. Spherical cap under appex load 
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Fig. 2. Load deflection curves 
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Fig. 4. Spherical cap subjected 
to uniform pressure 



Fig. 6. Cylindrical panel 
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Fig. 9. Deflection time history 
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applied load of 20 lbs/in is well below the linear critical (buckling) load for this geometry, which 
is 42.15 lbs/in. At a temperature of 50x F the shell is in a primary creep state for the first 28 
minutes, reaching a deflection of 0.2 in. and the critical time for creep buckling (this implies that 
the deflection becomes unbounded) is 35 minutes. At temperature of 500x F the shell ’maps’ into 
its post-buckled configuration almost immediately but the critical time for creep buckling remains 
almost unchanged. The dashed line in the figure represents a non-isothermal process where the 
temperature was suddenly increased from 50x F to 500x F after 0.3 hrs. As a result the shell 
snaps-through to its post buckled position at 500x F with small over-shoot and reaches its critical 
time of creep buckling two minutes sooner (33 min.). 

The last example is of shallow cylindrical panel made from Hastelloy-X and subjected to a 
concentrated force at the center of the panel (Fig. 8). Fig. 9 shows the deflection time history of 
points 1 and 2 to P = 0.6P c j which was kept constant. 

Discussion 

The study shows that in the presence of high temperature and viscoplasticity, the process 
of shell buckling assume an entirely new character. While the stability phenomena still exist 
under sufficiently large loads, buckling develops, as a time and temperature dependent deformation 
process under constant or variable thermomechanical loads of magnitude smaller than the purely 
elastic critical values. If the elastic behavior of a structure displays limit points and snap-through 
phenomena, the deformation process of creep buckling . become even more compli- cated and it 
usually exhibits a combination of snapping and creep responses. 
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